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ABSTRACT 


Fatogoma, Ouattara. MS., Purdue University, August 1990. Sampling Error Estimates 
of Satellite-Derived Tropical Rainfall Using General Circulation Model Data. Major 
Professor; Dr. Harshvardhan. 

The purpose of this work is to estimate sampling errors of area-time averaged 
rain rate due to temporal samplings by satellites. In particular, the proposed low 
inclination orbit satellite of the Tropical Rainfall Measuring Mission (TRMM; 35o 
inclination and 350 km altitude), one of the sunsynchronous polar orbiting satellites of 
NOAA series (98.890 inclination and 833 km altitude) and two simultaneous 
sunsynchronous polar orbiting satellites, assumed to carry a perfect passive microwave 
sensor for direct r ainfall measurements. This estimate is done by performing a smdy of 
the satellite orbits and the autocovariance function of the area-averaged rain rate time 
series. A model based on an exponential fit of the autocovariance function is used for 
acmal calculations. Varying visiting intervals and total coverage of averaging area on 
each visit by the satellites are taken into account in the model. 

The are generated by a General Circulation Model (GCM). The model has 
a diurnal cycle and parameterized convective processes. A special run of the GCM was 
at NASA/Goddard Space Flight Center in which the rainfall and precipitable 
water fields were retained globally for every hour of the run for the whole year. 

For the areas chosen( 5o by 4o grid boxes located at 130oE, 150OE, 160oW with 
latitude varying from 22oS to IQON), on average, the sampling error for the three 
months (December, January, and February) of the northern hemisphere winter would be 
of the order of 10% of the monthly mean for a satellite in the TRMM orbit. This error 
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is almost equivalent (difference less than 1%) to the sampling error produced by a 
polar-orbiting NOAA-series orbit. Observations with a system of two sunsynchronous 
polar-orbiting satellites simultaneously would reduce the sampling error from 10% to 
about 5%. 
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1. INTRODUCTION 

The roles of tropical rainfaU are very important in the energy budget of the 
planet and in affecting planetary waves that control weather patterns over the globe. A 
significant fraction of the heating of the tropical atmosphere comes from latent heat 
release in precipitating clouds. As a matter-of-fact, much of the sun s energy does not 
go into the air direcdy, but goes into heating land and ocean surfaces. Over the 
continents, the sun raises the temperature of the land surfaces which arc cooled as the 
air carries away the heat Over the ocean, particularly in the tropics, much of the sun's 
energy is used to evaporate water. It is only when this water vapor condenses to form 
cloud and rain that the sun's energy is deposited in the atmosphere. This energy, in 
turn, becomes the avadable energy that drives the winds. Since much of the tropics is 
covered by ocean, the formation of rain in the tropics is an important step in a process 
that transforms the energy in incoming solar radiation into kinetic energy which drives 
the motions of the atmosphere. Heated equatorial air moves towards the poles at high 
altirades and is replaced by cooler air flowing towards the equator at lower levels. 
Because of the vital role of latent heat released in driving general atmosphenc 
circulation, the variabitity of tropical rainfall is a key in understanding the circulation 

variability and short term climate changes. 

As we can see, tropical regions arc the primary source of the earth s weather and 

atmospheric circulation patterns. Accurate estimates of the amount of rainfall there 
would greatly improve our understanding of global weather, climate and the dynamics 
of the earth's atmosphere. However, tropical rainfall is not well monitored because of 
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its high variability in both space and time, and much of the area is covered by jungle 
and ocean. Consequendy, it is very difficult, even impossible to use conventional 
ground-based rainfall measuring devices such as rain gauges and radar. Because of this, 
rainfall observations in the tropics arc limited to what can be obtained from occasional 
ship and air reports, occasional experiments mounted specifically to collect data for a 
given region such as GATE [GARP Adantic Tropical Experiment, GARP (Global 
Atmospheric Research Program)], and extrapolations from land observations which are 
not easy to verify quantitatively. 

Continuous coverage of these large oceanic areas and highly variable rain rates 
is probably only feasible from space borne sensors which can provide a global data set 
over both the ocean and the continents. Measurements of cloud-top infrared emission 
from geostationary operational satellites carrying the Visible and Infrared Spin Scan 
Radiometer (VISSR) Atmospheric Sounder (VAS), currendy furnish the best indirect 
estimate of rain variability over the tropical region. In addition, passive microwave 
measurements made by the Electrically Scanning Microwave Radiometer (ESMR) 
aboard Nimbus-5, the Scanning Multichannel Microwave Radiometer (SMMR) aboard 
Nimbus-7, and the Special Sensor for Microwave/Imager (SSM/I) aboard U.S. Defense 
Department Satellites provide valuable information on rainfall rates and distributions. 
Some of these measurements however have been m ad e from geostationary satellites 
whose fields of view cover only a portion of the globe, and the estimating technique is 
based upon empirical procedure with fitted regression coefficients specific to the area 
and season of calibration. Besides, the estimates of average rain over a mouth have 
random errors greater than 50% of the mean. And the microwave measurements are 
maHe only from sunsynchronous polar-orbiting satellites. Due to the expected diurnal 
cycle of precipitation, measurements from sun-synchronous polar orbiting satellites 
would lead to a bias in the estimation of area-time averaged rainfall because they 
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revisit, at fixed local times, a given location. The low inclination orbits (LIO) or an 
observation system based upon two satellite orbits simultaneously would reduce 
significantly the observing time interval of an area down from twelve to six hours 
depending on the orbit parameters. Therefore, they could observe precipitation through 
the diurnal cycle. 

The use of satellites in measuring averaged rain rate within a given area during a 
period of time, introduces two significant sources of errors. Firstly, the error inherent in 
the method used to determine rain rate at a given instant. Secondly, there is the 
sampling error associated with the satellite only being able to observe a given area 
intermittently. The first source of error arises from a number of factors such as the non- 
homogeneous distribution of rain rates inside a field-of-view, the non-linear relationship 
between the rain rates and the radiant intensity (those two factors are known as the 
beam-filling problem) and calibration problems due to the limited information for some 
atmospheric and surface parameters. This error may in principle be considerably 
reduced by developing good retrieval algorithms to generate unbiased estimates of rain 
rate and designing sensor elements which help in calibrating satellite retrievals. 
Moreover, low-altitude satellites provide a s mall field of view which helps reduce errors 
due to the beam-filling problem introduced by the high spatial variability of 
precipitation. 

The second source of error, the sampling error, is due to temporal gaps in the 
observations induced by the orbit. Consequently, it is determined by the orbit of the 
satellite and the size of the swath scanned by the satellite as it passes over, and by the 
statistical characteristics of the observed rainfall. Its size depends on how well sampled 
the area is during a period of time. It is the dominant contribution to the error in time 
averages if the beam- filling problem and retrieval error are minimized. 
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The sampling error problems of tropical rainfall have been studied by many 
researchers mainly using data from cloud-top infrared temperature profiles and most 
importandy, the GATE rainfall data. The GATE experiment was conducted in two 
phases in the summer of 1974. The data were recorded from radars in an array of ships 
centered at 8.5°N, 23.5°W, off the west coast of Africa (Hudlow and Patterson, 1979). 
However, sampling errors may be different elsewhere on the globe or during a different 
season, since rainfall is known to be highly variable in both space and time. Laughlin 
(1981) examined the time series of area-averaged rain rates and noted that the lagged 
autocorrelation functions were close to exponential. His results show that for twelve 
hours repeated observations (Sun-synchronous Polar orbiting satellite) and for flush 
visits (satellite observes the entire box), the sampling error for monthly rain rate 
estimated over 2.5° by 2.5° grid box were about 8%. McConneU and North (1987) and 
Kedem et al. (1987) used a mixed probability distribution method along with an 
imaginary orbit ensemble to show sampling errors of about 8% and 10%, respectively. 
Shin and North (1988) used the random field method (Bell, 1987) along with uneven 
sampling intervals and fractional observations of 5° by 5° grid boxes to estimate the 

sampling error about 8% to 12%. 

We see that there are a number of sampling models. A survey of these sampling 
models tells us that most of them were developed not a long time ago and checked 
against the GATE rainfall data set only. The first of them is the poisson process. In 
this model, it is assumed that the individual field-of-view (FOV) of the satellite are rain 
fauge, the space-time statistics are homogeneous and stationary in the L Km by L Km 
averaging box and the averaging interval is 30 days. Moreover, the length and time 
scales are assumed to be short compared to these dimensions. An averaging box will 
have the results of N readings over the 30 days. If the probability of rain is P for 1 Km 
by 1 Km FOV and the average rain rate for pixels is R mm/hr and each FOV represents 
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a statistically independent observation (that is no FOV’s result of whether or nt it is 
raining depends upon the result found in a neighboring FOV) then we have a binomial 
process with N trials. One of the outcomes, rain is rare (P), hence we can use the 
poisson statistics approximation to the binomial distribution. The standard error of 
estimating the probability of rain expressed in percent is 1(X)/(P.N)°-^. This model 
underestimate the sampling error since the FOV's are not actually independent. 
Secondly, we have the Markov model of area average rain rate. This model, designed 
by Laughlin (1981) is based upon a stochastic model. The time series of area average 
rain rates is examined and if the lagged autocorrelation functions are reasonably close 
to exponential, we can constract a first order Markov process using the autocorrelation 
time corresponding to a particular averaging area (that is fit the autocoixelation 
functions to an exponential) and proceed to find an analytical expression for the 
sampling error. This is the model I used in this work against the GCM data. Thirdly, 
we have the partitioned rates and orbit ensembles model. This model was designed by 
McConnell and North (1987). First, an imaginary satellite orbit is started and flown 
over a given area at a time t, it returns 1 1 hours later and so on throughout the phase. 
During each flush visit, the individual pixels are collected into different categories 
(example: no rain, 0-5 mm/hr, 5-10 mm/hr, 10-20 mm/hr, and 20 mm/hr and above). 
The average rain contributed firom each of these categories is then estimated from the 
single imaginary orbit (realization). Realization number 2 is constructed the same way 
except that the orbit is started at time t plus one hour. In this way, an ensemble of 
satellite orbits can be constructed and the results can be compared. If the results differ 
widely from one orbit to another, it can be infeied for instance that significant storms 
are missed by the 11 hours gap. This model is so simple and assumption-free that it is 
very appealing. However, the autocorrelation in the area-wide time series of rain raes is 
not appreciated. Fourthly, we have the mixed distribution method. This model was 
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designed by Kedem, Chiu, and North (1987). It involves use of a sparse sampling in 
space and time on a regular grid. At the end of the phase, the pixel average rates are 
collected into a histogram of rain rates. These are then fitted to the log normal distribu- 
tion by a chi-square minimi zation procedure. After determining the appropriate 
probability density function (pdf) for the data, the average rain rate is calculated by 
finding the expectation value of the pdf. One such design of a satellite orbit 
corresponds to the sampling design of a satellite orbit withflush visits (exactly like 
McConnell and North above). Finally, we have the simulations with random fields 
model. It is designed by Bell (1987) following this recipe: first a gaussian random field 
is constructed infourrier spectral form over a grid. This can be thought of as a kind of 
turbulent vertical wind field. If the wind exceeds a certain threshold, we agree that it is 
raining. Such a procediue leads to islands of raining areas. The random variable 
representing the excedence over threshold can then be converted to a log normal distri- 
bution rain rate. T his model is more complete since the fields are constructed (it 
supplements the non-availability of real data) and evolved infourier space with conver- 
sions by way of fast fourrier transforms (a fourrier series can describe completely all the 
harmonics of diurnal cycles). Besides, Bell proved that once it is calibrated to GATE 
statistics, it imitates very well Laughlin's Markov model and the mixed distribution 
method. The data used in this snidy were generated using a General Circulation Model 
(GCM) that has a diurnal cycle and parameterized convective process (Randall et al., 
1989). The precipitation climatology of GCM is presented by Randall et al. (1989). 
Acmally, the data are obtained from a cloud generation process inwhich the vertical 
distribution of moisture content is parameterized. Throughout the year, roughly two- 
thirds of the simulated precipitation falls from cumulus ensembles, primarily in the 
tropics and the summer hemisphere; and the remaining third from large-scale saturation 
clouds, mainly in the winter midlatitudes. The geographical distribution of the 
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simulated precipitation is fairly realistic. The GCM successfuUy produces the 
rainbands associated with the Atlantic and Pacific branches of the intertropical conver- 
gence zone, and also the heavy precipitation of the south Pacific convergence zone and 
the equatorial western Pacific. The simulated seasonaUy varying precipitation maxima 
for eastern North America, tropical Africa, tropical South Amenca, and the Indian 
subcontinent are reaUsticaUy positioned, but somewhat strong. The model physics is 
updated every hour and usually only daily means are retained for analysis. A special 
run of the GCM was made at NASA/Goddard Space Flight Center in which the rainfall 
and precipitable water fields were retained globaUy for every hour of the run for the 

whole year. 

In the proposed work, I conducted a smdy of satelUte orbits along with the 
autocovariance function of the area-averaged rain rate time series. This was to estimate 
sampling errors for different grid boxes of 5° by 4° embedded in three different 
meridional regions; 150°E, 130°E, and 160°W (Figure 1). Note that the grid boxes at 
130°E and 160°W are completely in the ocean. At ISO^E we have two grid boxes; one 
in the ocean, the other over land. The orbit characteristics of three sateUite systems 
were used: the Tropical Rainfall Measuring Mission (TRMM) orbit, one NOAA-senes 
orbit, and the combination of two simultaneous NOAA-series orbits. 

TRMM is a joint Japan (Science and Technology Agency of Japan) / United 
States (NASA) mission which has set out to orbit satellites by 1994 designed to provide 
direct and quantitative measurements of rainfall data in the tropics (Simpson et al., 
1988). The main goal of the mission is to produce a monthly mean time senes of 
average rain rate over 5° by 5° boxes in the tropics. The TRMM orbit is planned to be 
circular with altimde 350 km and inclination 35° to the equator. This orbit gives exten- 
sive coverage in the tropics and allows extraction of the diurnal cycle in climatological 
rainfall because of its low inclination. The low altimde also ensures a small field of 



8 


view facilitating the interpretation of rainfall measurements over small spatial scales. 
Plans for TRMM include four principal sensors. The two microwave radiometers for 
r ainfall rates Over the oceans consist of a conically scanning multichannel instrument of 
moderate resolution with three dual polarized channels at 19, 37, and 85 GHz to provide 
a wider range of rainfall intensities, and a single-channel (19 GHz) electrically scanning 
instrument that provides a higher resolution at the more prevalent rain rates of 10 to 20 
mm hr*. The third sensor is a cross-track scanning precipitation radar, specially devel- 
oped for TRMM, which wiU measure rainfall over both land and sea. The radar data 
are necessary to obtain the vertical distribution of latent heat release, a very important 
parameter for global climate models. The fourth sensor is the Advanced Very High 
Resolution Radiometer (AVHRR) now in operation on polar-orbiting environmental 
satellites, which offers high-resolution channels in the visible and infrared spectral 
regions. The AVHRR will provide a link between measurements made by the first 
three TRMM sensors and those made simultaneously by visible and infrared radio- 
meters on geostationary spacecraft. 

The mission of the NOAA series is to collect global data on cloud cover, surface 
conditions such as ice and snow, surface and atmospheric temperature, and atmospheric 
humidity. It is also used to measure solar particle flux, and its third aim is to collect and 
relay information from fixed and moving data platforms. Finally, it provides continu- 
ous data broadcasts. An attempt is made to maintain two NOAA satellites in orbit at all 
times, a so-called "afternoon" bird with nominal observing times of 2 P.M. and 2 A.M. 
at altitude 833 km and "morning" bird with observations at 7 A.M. and 7 P.M. at 
altitude 870 km. The NOAA series is sun-synchronous with a circular orbit and 
inclination 98.89°. The satellite orbit of interest in this work is that of the latest 
"afternoon" bird NOAA- 10. The platform has four main sensors: the Advanced Very 
High Resolution Radiometer (AVHRR/2) composed of five bands; the Tiros 
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Operational Vertical Sounder (TOYS), a tiiree-sensor atmospheric sounding system 
composed of the High Resolution Infrared Radiation Sounder (HIRS/2) with a ground 
resolution of 17.4 km and 20 bands; the Stratospheric Sounding Unit (SSU) with a 
ground resolution of 147.3 km and three bands; the Microwave Sounding Unit (MSU) 
with a ground resolution of 105 km and four bands; the Space Environment Monitor 
(SEM) and the ARGOC Data Collection System (DCS). 

I also assumed that the satellite sees the whole grid box when it passes over and 
the rainfall is homogeneous inside a FOV. The orbit parameters and sampling intervals 
are obtained by a computer- generated model (Appendix) designed by Thomas Bell^ 
based on the theories outlined by D.R. Brooks (1977). 


1 Dr. Thomas L. Bell, NASA/Goddard Space Hight Center, Climate and Radiation 
Branch, Greenbelt, Maryland 2077 1 . 
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2. DERIVATION OF SAMPLING ERROR 


7 1 Gftneral Formulation 

As a measure of sampling error, I shall use the root mean squared (rms) 
difference between the actual area- averaged rain rate averaged over a month and the 
mt^an of the satellite observations during that month. The derivation of sampling error 
follows closely that outlined by Laughlin (1981). 

Let M(T) be the true monthly-averaged rain rate, 


M(T)= X 


T 
i ^ 


X(t)dt 


where X(t) is an area-averaged rain rate. 


( 2 - 1 ) 


X(t)=A (2-2) 

^A 

where R(r,t) is the rain rate at location r, t is the time since the beginning of the month 
and T and A are the period (1 month) and area averaged over. In practice, a continuous 
record of X(t) is seldom available. Assuming spatiaUy homogeneous rainfall statistics 
(so no weighting scheme is necessary), an estimate of M(T) by the satellite observation 
would be 

A N N 

M(T) = I fiXi(iAt) / S fj 
i=l i=l 


(2-3) 
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where N is the total number of observations. At the time difference between two 
observations (time difference between two consecutive nodes), T = N'At and Xj(iAt) 
would be a subarea averaged rain rate for the fraction of fi. 

As a first approximation, X(iAt) is substituted for Xj(iAt) because the rainfall is 
homogeneous in space and the satellite is expected to see the whole grid box. Then 
A N N 

M(T) = Z fiX(iAt)/ Z fj . (2-4) 

i=l i=l 


The mean squared enor (MSE) is the expected value of the squared difference of the 
two means. 

<ye2 = E[(M(T)-M(T))2]. (2-5) 

The sampling intervals are not always the same as the satellite rotates owing to orbital 
precession. Nonetheless, they would be different within only one nodal period. 
Moreover, ascending and descending nodes can be treated as two individual sampling 
series which have sampling time difference of At^ (visiting time difference between an 
ascending and descending node or vice versa). Then the estimated mean rain rate from 
satellite observation can be written as 


^ 1 
M(t)=rN N 1 


Z a; Z bj 

Li=l j=lJ 


[Z aiX(t-t-iAt) + Z bjX(t-KAtji -n jAt)] 
i=l j=l 


( 2 - 6 ) 


where aj and bj are the fractional coverages of each visit by the ascending and 
descending nodes, respectively, and t is time since the beginning of the month. 

I shall assume that if the satellite passes over a given grid box, it sees the whole 
box (flush visit). Therefore, a^ = bj = 1 and 
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M(t) = ~(2N) X(t+iAt) + 2 X(t+Atjj + jAt)] 
i=l j=l 


N 


The MSE <5^ = E[(M(T) - M(T))2] 


= E[M2(T)] + E[M2(T)] - 2E[M(T)M(T)] 


where 


E[M2(T)] =E[T2( X(t)dJ X(s)ds] 


(2-7) 


( 2 - 8 ) 


(2-9) 


Ao r_L_ N N 

E[MXT)] = EL(2N)2 { EX(t+iAt)} + { S X(t-i-AtH + jAt) } x 

i=l j=l 

N N 

{ L X(t-t-kAt) + Z X(t+At(4 -I- lAt) } ] (2- 10) 

k=l 1=1 


/s _1 N N . 

E[M(T)M(T)] = E[-(2N)-T] { I (t+iAt) + Z X(t-i-Au + JAt) } ’f X(t)dt] 

i=l j=l J 

o 


( 2 - 11 ) 


For a physical time -dependent function Q, its expected value E(Q) may be represented 
as the limi t of a time average of Q; if the limit exists. 

L 

( 2 - 12 ) 


E(Q) = ]im\^ I Q(a)da 


L— 


with L representing time. Therefore the true mean of X(t) denoted by E[X(t)] = X is 



13 


X = lim lT X(t)dt 
L->oc I 


The autocovariance function may be calculated as 
L 

R('C) = lim L ( [X(t)-X][X(t+x)-X]dt ^ 

L— 


where t is the time lag. Expanding the terms in the integrand yields 
L 

R(t) = lim l r [X(t)X(t+x)dt - ^ 

L->oo I 
-^o 

Let us evaluate the expected values of the expressions (2-9), (2-10), and (2-1 

L 

-ir JL 


L T T 


and 


E[M2(T)] = lim [q fx(t+a)dt fx (s+a)ds]da 

L-^ J J J 

O o o 

T T L 

E[M2(T)] = T^r r lim ["l r X(t+a)X(s+a)da]dtds 

J J L-^ J 


Letting u = t+a. 


T T 


E[m2(T)] 


f lim X(u)X(u-Hs-t)du]dtds^ 


L— >oo 


o *^o 


•^o 


(2-13) 

(2-14) 

(2-15) 

1 ): 

(2-16) 

(2-17) 

(2-18) 
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and 


= X2/ f [R(s-t) + X^Jdtds 


T T 


(2-19) 


if s-t = constant, so that R(s-t) = constant, we have a starionaiy time series. Therefore, 
we may simplify equation (2-19) through inspection of the area in the (s,t) plane over 
which the integration is performed. Choosing a coordinate system with one axis 
parallel to the lines s-t = c simplifies evaluation of the integral because R will be 
consant on all lines parallel to the axis. Let 
X = (s-t-t)/ 2, y = (s-t)/ JT 


Then, 


T/2 x = 2T-y 

E[M^(T)]=:^ r r [R (^y)-hX2]dxdy 
■^o -^x = - 2T-i-y 


Letting U =j2y, 


E[M2 (T)] =4/ ( 1 -UA’)[R(U)-t-x2 ]dU 


( 2 - 20 ) 


( 2 - 21 ) 


E[(M2(T)] = lim l r ^2 ^ X(L^-iAt)X(L^-lcAt)]dL 


L — >00 


i=l k=l 


L 

+ ^“°L (2N)2 [I 2X(L-HiAt)X(L-hAtd-HlAt)]dL 

L — >00 j i=l 1=1 

^o 

L 

if _L_ N N 

L (2N)2 [Z Z X(L-HAtd-HjAt)X(L-i-kAt)]dL 

L->oo / j=l k=l 


o 


15 


L 

1 T-L- N N 

+ lim L (2N)2 [Z I X(L+Atd+jAt)X(L+Atd+lAt)]dL 
j=l 1=1 

( 2 - 22 ) 


Working out the integrals in (2-22) by almost the same procedure as above. 


E[]^^(T)] = (2N)2 


N-1 

2NxR(0)-i-2NxR(At_,) -t- 2 S 2(N-k)R(kAt) 
^ K=1 


N-1 

+ 21. (N-k)R(kAt-i-Atd) 
k=l 


N-1 

-H 2 S (N-k)R(kAt-Atd) 
k=l 


(2-23) 


L 

E[M(T)M(T)] = iiin 


N 


(2N)T [Z X(L+iAt) X(L-i-t)dt]dL 

i=l J 

o -^o 


+ lim 
L— 


'' o 


N T 

(2N)T [Z X(L+At(i-HjAt) r X(L-Ht)d]dL 


(2-24) 


Using the same procedure, 


E[M(T)M(T)] = (2N)T 


T 

r _N 


N 




[z [R(t-iAt)-hX2]dtH- Z [R(t-Atd-jAt)-^X2]dt] 

i=l j=l J 

o ^o (2-25) 


Finally, from (2-8), (2-21), (2-23) and (2-25), 
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T 


o| = ^ r (1- -f) R(x)dT 

Jo 

+ (2N)2 { 2NR(0) + 2NR(Atji) 


N-1 

+ 22 2(N-k)R(kAt) 

K=1 

N-1 

+ 2 2 (N-k)R(kAt+Atd) 
K=1 

N-1 

+ 22 (N-k)R(kAt-Atd) } 
K=1 


- 2 _ 

■ (2N)T 


N 


{I 

k=l 


T 

r 



RCt-kAOdT + ^ f R(x-kAt-At(j)dT ] 


(2-26) 


2.2 Sun-sv nchronous Polar Orbit 

For a sun-synchronous polar orbit, inserting At^j = At/2 = 12 hours in the general 
expression (2-26) suffices to describe the sampling error. 

2.3 The Case of Two Simultaneous Sun-sv nchronous Satellite Observations 

In observations by two simultaneous SSO satellites, it is assumed that 
precipitation is observed by the same microwave sensors and they are perfect so that no 
optimal weighting scheme is required due to different measurement errors for the 
different spacecrafts. 

The two satellite observations can reduce the observation interval to shorter than 
12 hours depending on their phase difference Atj. If At^j is six hours, four observations 


17 


are available in a day for a given area with 6 hour intervals. In this case, the estimated 

mean can be expressed by: 

/\ JL r N N , 

= (2N) L ^ X(t+ito) + ^ 

i=l j=l ’ (2-27) 

where tg is 12 hours. 

Following the same steps of calculations as in chapter 2.1, the MSE for 2 SSO satellite 
observations is the same as (2-26) except At is replaced by Iq, and their difference is in 
the values of At, Atj, and tg. 

I shall express the sampling error by standard error, that is in percent of monthly 
area-averaged mean rain rate (lOO’Og /X). 
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3. SAMPLING CHARACTERISTICS OF A SATELLITE 

3. 1 Orbital Dynamics for Sampling Error Study 
The sampling error defined in equation (2-26) depends on the sequence of 
satellite observing times and sampling intervals, since this sequence determines the 
accuracy of the estimated mean rain rate from satellite observations. In generating 
sequences of satellite overpasses, it is assumed that the orbits are circular. 

According to Brooks (1977), the longimde-latimde history relative to an initial 
point on the equator is calculated as 

(}) = sin'^[sin (i) sin((0-i-0)t)] (3-1) 

X = tan-l[cos (i) tan((0+0)t)] - (O-co)t (3-2) 

where 

(j>, X latimde and longimde, 
i inclination of satellite orbit, 

0 perturbed mean angular velocity of the satellite on its orbital plane, 

0 angular velocity of precession on its orbital plane due to earth’s oblateness, 

0) angular velocity of precession on the equatorial plane due to earth's 
oblateness, 

Q earth's angular velocity on the equatorial plane, 
t time, and 
h the satellite altitude. 

Let 0 + 0 = a and - co = p. Then from (3-1) and (3-2), we have 
t = ^ [2:tn + siriK^t^] = nXj^ + 1/a siri^ 


(3-3) 
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X, = -§ [2m + sin*^(^*i)] + tan ^cos i tan(sinH^|»]_ 


(3-4) 


Where n (integer) is the number of rotations and Tn is the nodal period of the satellite 
for the given orbit. Then (3-3) gives the time at which the satellite crosses latitude (|) in 
rotation n, and (3-4) gives the longitude when the satellite passes latitude 0. 

Xjj = 2 jc/o. 

For the ascending node, 



(3-5) 


while for the descending node. 


td = a[2ra + 7t-sin-iC^‘*i]. 


(3-6) 


Sinulariy, 

= -|[27cn + siri^ (^f)] + tan ^[cos i tan(siri^ (S^^l and (3-7) 

= -|[27in + 3t - sin ^ (s^^] + 7t - tan ‘Hcos i tan(sm^ . (3-8) 


From (3-5), (3-6), (3-7), and (3-8), 

h ■ ^a ^d , ' ^d = iTz/a = '*'n, (3-9) 

n+1 n n+1 n 

- ^a = ^d ■ ^d = = -pxn^ (3_ 10) 

n+1 n n+1 n 


and 



,X = p/a [tz - 2sin^ + ^'^tan 

n % 



(3-11) 

(3-12) 
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Therefore, the ascending and descending nodes can be treated as separate sequences 
which differ in initial phase depending on <j). Define the repetition factor 
Q = ©+6/S2-(o = ct/p 

Thus, 

. f ^aj} * Q = 2:t/p = f(i,h) = constant, for a given i,h, and 

(3-13) 

^^d ~ (^d^' ’ Q = [7t-2sin (sin<|)/sin i)] = f(i,h,0) = 

constant, for a given i, h, 0, (3-14) 

At is the time from one observation by an ascending/descending node to the next 
observation by an ascending/descending node over a particular spot on the earth. 
Likewise, At^j is the time from one observation by an ascending/descending node. For 
the sun-synchronous oihit. At is equal to one day. If the orbital parameters do not 
satisfy the aforementioned condition, that is for a low inclination orbit, the orbital plane 
will process either eastward (At > 1 day) or westward (At < 1 day) with respect to the 
earth. The deviations of At from one day for various orbits are summarized in figure 2. 
In this figure, the difference is in minutes between the visiting interval of the given low 
inclination orbit and that of a sun-synchronous. The visiting interval by an 
ascending/descending node wdl be limited to within 30 minutes of 24 hours. 

12 Sampling Characteristics o f the Satellites 
I shall assume that the viewing characteristics of the satellites are similar to 
those of the Electrical Scanning Microwave Radiometer (ESMR) flown on the Nimbus- 
5 satellite. Even though this sensor is not aboard NOAA, I shall assume that it is. In 
fact, I consider the orbit of NOAA- 10 with an imaginary ESMR (the one planned to fly 
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on TRMM satellites) flown on it The ESMR measures the intensity of microwave 
radiation at 19.35 GHz in 78 fields of view (FOV's) during each scan of the instrument. 
The instrument scans perpendicularly to the direction of satellite motion, from about 
50® to the left of nadir to about 50° to the right of nadir, requiring 4 seconds per scan. 
The beam widths of the FOV's ranged from 1.4° x 1.4° near nadir to 1.4° x 2.2° at the 
50° extremes. This sensor is interesting in this study because its FOV's overlap 
substantially, and nearly complete coverage of the area within the view of the satellite is 
assured. This condition is necessary for the assumption of flush visits. Consequently, 
for a better and total coverage of the 5° x 4° grid boxes, I extend the average effective 
size of the swath to about 54° on either side of nadir. 

Knowing the scan angle (54°), the longitude and latitude of the grid box center, 
the altitudes of the satellites, and the size of the grid box (5° by 4°), the parameters of 
the orbits and the visiting times of the satellites are calculated using the computer model 
in Appendix. Table 1 summarizes those parameters for the TRMM orbit and the 
NOAA orbit satellites assumed to carry the ESMR. Figure 3 presents the variations of 
sampling interval At^j versus latitude for the TRMM orbit This variation is caused by 
the orbit precession. For the sun-synchronous orbit, there is no orbital precession and 
At and At^j equal 12.00 and 6.00 hours respectively. Recall, for given inclination i, and 
altitude h. At is constant while Atj is function of latitude (j). 
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4 . RAINFALL DISTEUBUnON ESr TIME 

4.1 Examples of Rai nfall Distriburion 

Figure 4 gives the distribution of area-averaged rain rate over the grid boxes 
considered during the months of December, January and February for the GCM data. 
The monthly mean rain rates and their standard deviations are tabulated in Table 2 and 
presented in Figures 5 and 6. As can be seen, area-averaged rain rates are highly 
variable with respect to time and gnd box location. Although the period considered 
corresponds to the rainy season of the region (south of the equator), the monthly 
averages of area-averaged rain rate are moderate. They vary from 0.08 mm/hr to 1.13 
mm/hr during the three months. A close look at the rain rate distribution shows the 
existence of diurnal cycle, especially the rain rate distribution in the grid box centered 
at 150°E, 10°S. It should however be indicated that I did not smdy the temporal 
variability of rain rate. 

4.2 The Autocovariance F unctions 

The autocovariance function is determined by investigating the GCM area- 
averaged rain rate time series and using the formula 
L 

R(t) = lim L r [X(t)-jq[X(t+T)-jqdt 
L-^ I 

(4-1) 

The autocovariance functions calculated for different grid boxes for the three 
months of December, January and February are presented in Figures 7, 8, and 9. All the 
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derived curves almost appear to be examples of exponential decay. Each is least 
squared fitted to a function of the form 
,x . 


R(x) = B exp(- 1 Xq I ) (4-2) 

where B is the coefficient of variation and Xq is the autocorrelation time. This fit is 
simply a first-order Markov process. The best fit values of B and Xq are tabulated in 
Table 2 and Figure 8 presents the variation of the correlation time versus grid box 
center. 


4.3 Final No te on the Sampling Error Derivation 
The analytic form (4-2) of the autocovariance function allows for simpler 
evaluation of the sampling error using the following relationships: 

T -T/t„ 

r RCDdT = BXo(l-e ) (4-3) 

> o 


xR(x)dx = B(Xq2-x 2 e 


-T/x, 


-T/x, 


-Tx^e 


o 


(4-4) 


R(x-kAt)dx = 2 BXo-BXq2 


-kAT/Xg -T/Xq kAt/Xo 
e - BXqO e 


(4-5) 



R(x-kAt-At(j)dx = 2 BXq-BXq2 


-T/Xq kAt/Xo At(j/Xo -kAt/Xo At^j/x^ 

e e e - e e 

(4-6) 


and the geometric series 


N i Afl-A^i 
2 A = (1-A) 
i=l 


(4-7) 
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Using (4-3), (4-4), (4-5), and (4-6) in (2-26) yields 


At^/To Ai/Xq -Atji/Xo 

2 -J r T f (e +l) + e lie 

= (T/Xo) + 2x„ (gAt/Xo . 1) 


±H) 


-T/Xq Atjj/Xo At/Xo -Atyxo 

At {£s + 1 ) + g (s 


+JQ) 


At/Xo -Atyxo -At(i/Xo 

A ts t., 2 ±S 

+ 2 V (^At/T„ . 1)2 



(4-8) 


For observations with one sun-synchronous polar orbiting satellite, At^j^ = At/2 = 12 
hours. With simultaneous measurements from sun-synchronous polar-orbiting 
satellites. At = to = 12 hours and At^^ = 6 hours. 
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5. SAMPLING RESULTS AND DISCUSSIONS 

Table 2 summarizes the sampling errors of the monthly area averaged rain rate 
estimates for the three different orbits of TRMM, NOAA and the 2 simultaneous sun- 
synchronous orbit over the three months (December, January, and February) of the 
Northern Hemisphere winter. Figures 9 to 20 present the variation of sampling errors 
of the monthly area averaged rain rate versus gnd box center, monthly mean area 
average rain rate, standard deviation of the rainfall and correlation time of the rainfall 
time series. 

In December, on the average, the sampling error for the monthly mean rain rate 
estimate over a 5° x 4° grid box is 9.33% for the proposed TRMM orbit satellite 
observations, 9.61% for the NOAA orbit sateUite observations and 4.33% for 
simultaneous observations with two sun-synchronous polar orbiting satellites. 

In January, it is 9.94% for the TRMM orbit satellite observations, 10.25% for 
the NOAA orbit satellite observations and 4.86% for simultaneous observations with 
two sun-synchronous polar orbiting satellites. 

In February, it is 10.15% for the TRMM orbit satellite observations, 10.46% for 
the NOAA orbit satellite observations and 4.95% for simultaneous observations with 
two sun-synchronous polar orbiting satellites (see Figure 21). 

As an observation, sampling errors for the three orbits increase somewhat from 
December to February. This is mainly due to the correlation of rainfall; the average 
correlation time decreases substantially from about 13.30 hours in December to about 
10.00 hours in February. At the same time, the monthly mean, the standard deviation of 
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the area averaged rain rate and the coefficient of variation of the autocovariance 
function remain almost constant from December to January, and decrease from January 
to February (see Figure 22). The forms of the variations of sampling errors versus 
monthly rainfall mean or autocorrelation time or standard deviation of the rainfall 
distribution for the three orbits during the three months are similar. The differences, on 
the other hand, are in the magnitudes of the different parameters indicating that the 
overall area-averaged rain rate statistics are similar for given seasons and areas. As a 
consequence, one may use the same statistical model to study rainfall in an area during 
a season. Moreover, sampling errors have a tendency to decrease with increasing 
monthly mean area-averaged rain rate and to increase with increasing standard 
deviation because the less rainfall the satellite senses, the less accurate is the 
measurement. The more variable the rainfall distribution is in time, the less accurate is 
the measurement by the sensor. The non-linear decrease and pseudo wave--Uke 
variations depend on the contribution of the autocorrelation time Tq in the formulation 
of sampling error (see formula 4-8). The value of this parameter is not constant, 
therefore, its combination with exponentials affects greatly the sampling errors giving 
them sometimes pseudo-periodic variations. The visiting intervals At, is constant for a 
given orbit characteristic (inclination and altitude). The sampling intervals Atjj on the 
other hand decreases linearly with increasing latitude. The wave-like variation of the 
sampling error against the autocorrelation time and latitude results from the same 
reason. 

The sampling error of the TRMM orbit is almost equal to that of the SSO 
NOAA. The satellite orbit study shows that for grid boxes in the equatorial region, as it 
is in the case of the present study, the visiting time interval of a low inclination orbit 
satellite is almost equal to that of the sun-synchronous orbit. In Fig\u-e 3 we see that the 
sampling interval varies from 11.68 hours to 11.83 hours for the range of latitude 
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between 22°S and 14°N. This is not very different from the 12 hours of the sun- 
synchronous orbit. Typically, there are no more than two observations a day. As a 
result, the sampling errors are almost equaL Away from the equator, a low-inclination 
orbit satellite may revisit a given location more than twice a day. It follows then that 
the visiting interval reduces and so does the sampling error compared to that of the sun- 
synchronous orbit. Even though those two orbits provide almost the same sampling 
error, the FOV of the TRMM is smaller than that of the sun-synchronous orbit because 
of its low altitude (in the present study almost half of that of the sun-synchronous). 
Consequendy, the error due to the beam-filling problem would be reduced. 
Furthermore, the expected additional errors due to the diurnal cycle of the rain field will 
be added to the sun-synchronous orbit because of its fixed visiting sequences, twice a 
day, while they tend to be cancelled by the TRMM orbit Because of this, the overall 
error would be smaller with a TRMM orbit. 

The sampling error reduces by a factor almost of two when considering 
simultaneous observations with two satellites sun-synchronous orbit. In the same vein, 
the diurnal cycle in the rain field would not be significant in this case because with four 
observations a day, the system can resolve all the harmonics of the diurnal cycle. Yet, 
the large field of view, due to the high altimde, could add large errors caused by the 
beam-filling problem. 

When the monthly mean area-averaged rain rate is high, for instance, greater 
than 0.4 mm/hr, and when the autocorrelation time is less than or equal to 12 hours, the 
sampling error is roughly about 8% for the TRMM and sun-synchronous orbits and 3% 
for the two simultaneous sun- synchronous. This result agrees with the former studies 
done Avith the GATE data assuming flush visits and varying visiting time intervals (Shin 
and North, 1988). 
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6. SUMMARY AND RECOMMENDATIONS 

A study of satellite orbits has been conducted to estimate sampling eixors of 
area-time averaged rain rate caused by gaps in sateUite observations. This study 
focused on the proposed low-inclination orbit of TRMM, the sun-synchronous polar 
orbit of NOAA-10 and simultaneous measurements from two sun-synchronous orbits. 
This estimate was done by performing a study of the autocovariance functions of the 
area-averaged rain rate time series. A model based on the exponential fit of the 
autocovariance function is also used for actual calculations. Varying visiting intervals 
and total coverage of the averaging area on each visit by the satellites were taken into 
account in the model. 

Of the areas chosen (5° by 4° grid boxes located at 130®E, 150°E, 160°W with 
latitude varying from 22°S to 10°N), on average, the sampling error for the three 
months (December, January, and February) of the Northern Hemisphere winter, would 
be about 10% for a TRMM orbit. This error is within 1% of the sampling error 
produced by a polar-orbiting NOAA orbit The systematic error is higher in sun- 
synchronous polar-orbiting sateUites because of their large fields of view. On the other 
hand, sateUites in low-inclination orbits have the potential of filtering expected diurnal 
cycles from the observation because of their relatively low orbit and inclination. Thus, 
the preceding findings indicates that tropical precipitation observations with low 
inclination orbits are more efficient. A simultaneous observation with a system of two 
sun- synchronous polar orbiting satellites would reduce the sampling error from 10% to 
about 5%. Such a system would be very effective in observing precipitation which is 
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very variable in space and time insofar as sampling error is concerned because its four 
or more observations a day can average out all the harmonics of the diurnal cycle. Still, 
the large field of view may introduce a large systematic error. 

The results obtained in this work provide a quantitative and pedagogic basis for 
estimating sampling errors as far as the GCM data is concerned. The global results will 
apply only where the time scales (particularly the diurnal cycle) of the precipitation and 
the fractional visiting areas are simulated realistically by the model. Another question 
seldom touched on so far is the idea of merging data from different spacecrafts. Each 
will have its own measurement error structure as well as its own sampling 
characteristics. An optimal weighting scheme will have to be devised for the 
integrations. Microwave data from polar platform can be merged with the low 
inclination orbiter data by some optimal procedure. We also need to understand the 
contribution that can be made by integrating other satellites such as GOES. The data 
used in this work are simulated data. Also, the previous smdies depended upon a single 
oceanic tropical rainfall data set. As a consequence, we must be concerned about their 
representativeness for conclusions. For rainfall statistics characteristic of the Inter 
Tropical Convergence Zone (ITCZ) conditions, the results arc in line with those 
obtained with the GATE data for the same assumptions. However, one should be 
concerned about the uncertainties of some parameters in the statistical model such as 
the values of the autocorrelation times and coefficients of variation extrapolated from a 
least-square fit process. 



Table 1. Orbits parameters for TRMM and NOAA-10. 


ORBIT 

TRMM 

NOAA-10 

ALTITUDE 

(kilometers) 

350.0 

833.0 

INCLINATION 

(Degrees) 

35.0 

98.9 

SCAN ANGLE 
(Degrees) 

54.0 

54.0 

ALPHA 

(Radian/Second) 

0.00114681 

0.00102983 

BETA 

(Radian/Second) 

0.00007429 

0.00007272 

NODAL PERIOD 
(Minutes) 

91.31 

101.69 

PRECESSION OF LONGITUDE 
(Degree/Orbit) 

-23.32051 

-25.04531 

SWATH HALF-WIDTH 
(Kilometers) 

510.37 

1353.63 

DELTAT AT 
(Hours) 

23.49 

24.00 



Table 2. Rainfall statistics for December, January, and February. 
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Table 2, continued 
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Table 2, continued 
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Figure 1. Grid box map. 
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Figure 4. Distribution of the monthly area averaged rain rates versus time for different grid box locations in December, January 
and February. 
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Figure 22. Variation of averaged rain rate statistics from December to February. 
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#include<stdio . h> 
#lnclude<math . h> 
♦include " out_hdr . h " 
♦include "orb_func . h" 
♦include "orb_func . c" 
♦include "orb_plot . c" 


void main ( ) 

/ 


/•k^-k-k-k-kiritiicirir-kie* 






int whichOne; 
BOXOBS */ 


/* Set WhichOne = NUMOBS, WHERESAT, or 




*★★*★★*★* 


if ( setupOutput ( ) ) goto error; 


WhichOne - BOXOBS; 

switch ( WhichOne ) { 
case NUMOBS: 


{ 


if { NumObservO ) ( 

P^iiitf ("\nProblem in NumObserv ! " ) ; 
goto error; 

} 

break; 


} 

case WHERESAT: 

{ 

WhereSat ( ) ; 
break; 

} 

case BOXOBS: 


{ 


if ( BoxObservedO ) { 

f {"\nProblem in BoxObserved ' " ) • 
goto error; 

} 

break; 


error ; 


GLOBAL VARIABLES 


*/ 

int direction, printLevel=l, boxMatters, retrograde, symmetryUsed; 
double 

Pi, TwoPi,Pi2,DtoR, RtoD, incR, sininc, cosine, alpha, beta, betaOalpha; 

double maxLatTr, delta, sinDelta, cosDelta, retroPi, tauN; 

double 

t, longSat, travelTime, crossLong, crossTime, dt, tStart, tLim, longBoxR; 

double phiN,phiS, longW, longE, tl; 

double longLeft, windowLong, dTdL [2] , tLeft [2] ; 

/* 'printLevel' indicates how much printing is allowed. 

'boxMatters' indicates how much box size needs to be specified, 
'retrograde' is set to 1 if a retrograde orbit, i.e., inc >= 

90! . 

'symmetryUsed' is set to 1 if box center is south of Equator. 

*/ 

/★★*★******* 

*★★★★★**/ 


int setConstants ( ) 


/* Sets up constants needed for orbital calculations 

and prints out some interesting 


orbital parameters */ 


{ 

int i,direct0, swathWidthInKms,boxSizeInDegs; 

double alt, inc, scan^phiBox, longBox^delx, dely, swathWidth; 

double rEarth, ThrHf Jo2^mu; 

double 

rSat^ reOrsat, tau0,M0,Mdot, omegaDot, omegaEarth, omegaCapDot ; 
double longShift, temp, phiSatO, longSatO; 
double phiBoxR, phiTrack, longTrack, side, delLat , delLong; 

/* 

★ Initial conditions 


* 


*/ 


it 


Define parameters of orbit and scan angle 
alt =870.0; /* altitude of satellite in kms */ 

inc - 98.89; /* inclination of orbit plane in degrees 

*/ 

swathWidthInKms - NO; /* Next line ignored if NO */ 
swathWidth - 1200.; /* Side-to-side swath width in kms 

scan = 54.0; /* Scan angle in degrees. Ignored if 

swathWidth- 
InKms - YES. */ 

/* Define location of box observed */ 
phiBox - 14.0; /* Latitude of box center in degrees */ 

longBox = 130.0; /* Longitude of box center */ 

/★ Define size of box. */ 
boxSizeInDegs = YES; 
delx = 512 .; 


/* width of box, kms */ 
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dely = 512.; I* height of box, kms */ 

I* The next two specifications are used if 

boxSizeInDegs = YES. */ , . 

delLong - 5.; /* width of box in degs. */ 

delLat - 4.; /* height of box in degs. */ 




I* Define initial conditions for satellite */ 
phiSatO « 0.; 


longSatO * 0.; 
t - 0.; 

direct 0 - 1; 

descending orbit, 

ascending orbit */ 
tLim « (30,) * 86400,; 
orbit */ 

dt • 1 . * 60 . ; 
posit, every dt (secs) */ 
tStart “ 0.; 
printing satellite 


/* initial latitude 
/★ « longitude */ 


/★ ” direction; 

direction = 0 for 

direction - 1 for 

/* Time (secs) allowed satellite to 

/* WhereSatO shows sat. 

/* Time when WhereSatO starts 
position. */ 


if( printLevel > 0) { ^ , v „ 

SHOW "Satellite altitude = %6.1f kins.Xn ,alt); 
SHOW ” inclination » %6.1f 

degrees . \n”, inc) ; r, x 

if ( ! swathWidthlnKms ) 

SHOW "Scan angle - %6.1f 

degrees . \n\n” , scan) ; 

if ( boxMatters ) { 

if ( boxMatters ~ BOXOBS ) { 

SHOW "Latitude of box center * %6.1f 

degrees. \n",phiBox) ; . . ^ w 

^ SHOW "Longitude of box center » 

%6.1f degrees. \n", longBox) ; 

if ( boxSizeInDegs ) { 

if ( boxMatters — = BOXOBS ) 

SHOW "Box width » %6.1f 


degs . \n" , delLong) ; 
degs . \n\n" , delLat) ; 

]cms.\n",delx) ; 
kms . \n\n" , dely) ; 

} 


SHOW "Box height = %6.1f 


} 

else { 


if ( boxMatters BOXOBS ) 

SHOW "Box width » %6.1f 

SHOW "Box height * %6.1f 


} 


if ( boxMatters !* NUMOBS ) { 

SHOW "Initial conditions for 


satellite :\n") ; 


degrees . \n" , phiSatO) ; 
degrees . \n" , longSatO) ; 


SHOW "\tLatitude - %6.1f 
SHOW "XtLongitude * %6.1f 
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SHOW "NtTime - %6.1f secs.\n",t); 

SHOW "\tdlrection = %d. \n\n",directO) ; 
SHOW "\tT±me Limit - %6.2f 

days . \n\n " , tLim/ 86400.) ; 


constants 


*/ 

*/ 


/*Now do a lot of calculations and introduce physical 
necessary for calculations */ 


Pi - PI; 

TwoPi “ Pi + Pi; 
Pi2 = PI2; 

DtoR = Pi/180.; 

RtoD - 180. /Pi; 


/* Pi / 2 */ 

/* for converting degrees to radians 
/* ” ” radians to degrees 
/* radius of earth, )cms */ 


/* Three halves J/2, e's quad 


rEarth - 6378.145; 

ThrHfJo2 “ 0.0016238235; 

moment */ */ 

mu ■■ 398601.2; /* gravitational const of earth / 


rSat « alt + rEarth;/* radius of sat. orbit */ 
incR « inc * DtoR; 
sininc - sin(incR); 
cosine - cos(incR); 

retrograde * (cosine <* 0. ) ? YES : NO; 
maxLatTr - ineR; ^ 

if ( retrograde ) maxLatTr * Pi * incR; / Maximum orbit 
1 at it * / 

’retroPi - ( retrograde ? -Pi : Pi ) ; /* Used in getting 

longitude — 0. ) retroPi - 0.; /* from alphat. */ 


reOrsat * rEarth / rSat; 
if ( swathWidthInKms ) 

delta * .5 * swathWidth / rEarth; 

else 

^ temp =* (rSat/rEarth) *sin (scan*DtoR) ; 

if (temp >1.) temp - 1.; /* check if view extends 

beyond earth*/ 

delta - asin(temp) - scan*DtoR; 

^ /* swath width in radians */ 

sinDelta - sin (delta); 
cosDelta - cos (delta); 


tauO - TwoPi * rSat * sqrt(rSat / mu) ; 


/* unperturbed 
Z’* unperturbed mean motion, 


period, secs */ 

MO - TwoPi / tauO; 

rads /sec * / , 

Mdot = MO * (l.+ThrHfJo2 * reOrsat*reOrsat* (l.-1.5*sinlnc * 

sininc) ) , perturbed mean 

motion, rads/sec */ 

theta * / 

omegaDot - ThrHfJo2 * reOrsat*reOrsat * Mdot * (2.- 

2 .5*sinlnc*sinlnc) ; 
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/* perigee 

velocity# rads/sec */ -i • ^ */ 

tauN - TwoPi / (Mdot + omegaDot) ; /* nodal period^ secs / 

omegaEarth - (360.9856473 / 86400 .) *DtoR; /* earth rotation 

rate^ 


rads /sec ^ / 

omegaCapDot - -ThrHfJo2 * reOrsat*reOrsat* ^cosine 

orbital plane, rad/s*/ 

alpha - Mdot omegaDot; 

beta - omegaEarth - omegaCapDot; 

betaOalpha « beta / alpha; ^ 

longShift * - beta * tauN; /^change in satellite long. 

after 1 orbit*/ 

crossLong * longShift; 
crossTime * tauN; 


if ( printLevel > 0) { . « ^ 

SHOW "Swath half-width * %5.2f degreesXn , delta RtoD) ; 
SHOW "Swath half-width - %5.2f kms\n\n",delta*rEarth) ; 
SHOW "Nodal period = %9.4f mins . \n", tauN/ 60 .) ; 

SHOW "alpha - %10.8f radians/sec. \n", alpha) ; 

SHOW "beta » %10.8f radians/sec. \n", beta) ; 

SHOW "beta/alpha * %l0.5f .\n", betaOalpha) ; 

SHOW "precession of longitude =* %10.5f 

degrees/orbit . \n\n" , longShift * RtoD) ; 

) 

if ( ! boxMatters ) , ^ j 

goto endorbit; /* Following is needed for running 

BoxObserved 

and NumOb s e r v * / 

/* BOX MATTERS */ 

/* Determine boundaries of box observed */ 

if ( boxMatters == NUMOBS) { /* If NumObserv is being run 

*/ 

phiBox - 0 . ; 
longBox * 0 . ; 
delLong - delx * 0 . ; 

) 

phiBoxR - phiBox * DtoR; 
longBoxR - longBox * DtoR; 
if { boxSizeInDegs ) 

temp * 0.5 * delLat * DtoR; 
phiN * phiBoxR + temp; 
phis - phiBoxR - temp; 
temp - 0.5 * delLong * DtoR; 

else /* box size specified in kms */ 

temp » 0.5 * dely / rEarth; 
phiN » phiBoxR + temp; 
phis = phiBoxR - temp; 

temp = 0.5 * delx / (rEarth * cos (phiBoxR) ) ; 
if ( boxMatters =* BOXOBS && temp > Pi ) { 
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printf("\n The box is too wide."); 
return ( 1 ) ; 


) 


longE “ temp; 
longW « -temp; 

symmetryUsed * NO; 

if ( phiN < fabs(phiS) ) 

requires */ 

{ 


/* Check whether box position 
/* reflection 


symmetry to be used. */ 
temp ” -phiN; 
phiN - -phiS; 

temp; t-his case all latitudes must be */ 

directO - -dxrectO; /* In this case axx */ 

' /* output in BoxOPsurv 

with opposite signs. */ 

) 

switch ( i • boxCheckO ) { 

case 1 : { , ,, v . 

printf("\n Box height too large. ), 

return (i) ; 

} 

printf("\n Box outside latitudes viewed by 

satellite."); 

return (i) ; 


1 


/* Box seen every orbit! */ 


case 3 : 

^ printf<"\n Box is at least partially viewed 

every orbit."); pj.j_ntf("\n Unable to con?>ute portion seen each 

time."); 

return (i) ; 


1 


} 


/* Since the box long, has been moved to 0, the 
satellite^has to^ also shifted. Note that in outputting 

longitude, longBox should be added. */ 

longSatO — longBox; /* Note that this does not affect 
WhereSat */ 

if { /^J^atellite position doesn't matter to 

NumObserv */ 


/★ F.ND OF BOX MATTERS 


★ / 


endorbit : 


/* follow satel. to Eq. and get its longitude 


there ...■*/ 

direction * directO; 

getLongO ( phiSatO*DtoR, longSatO^DtoR, 1- 
direction^ 0 . , filongSat ) ; 


been moved 


/* Note that if satellite was descending, it has now 
to longitude at equator where it is ascending. 


lonaSat * Normalong { longSat ); / 
t +- travelTime; /* Set time to time of equator crossing V 
direction - 1; /* Satellite is at equator, ascending, 

orbiting. */ 


return (0) ; 


★ 

* 


★ 

*/ 


int boxCheckO 

{ /* Checks how box is observed by satellite and reports 

result */ 

if(phiN>Pi2 II phiS<-Pi2 ) 

return (1) ; 

/* Now check whether box is never observed, or always. 

^ if ( phis > maxLatTr + delta ) /* never seen. */ 

return (2) ; 


1 1 


delta) ) 


if ( (maxLatTr + delta < Pi && phiS <- - maxLatTr + delta ) 
(maxLatTr + delta > Pi && phiN >=■ Pi - maxLatTr 


return (3) ; 


/* Box seen every 


orbit */ 

return (0) ; 

} 


/* 

★ 


* 

*/ 


void orbitSat ( drawPtrO, frcPtrO, longlist,nCross,nWindows,draw ) 

DRAWTEMPLATE *drawPtrO; 

FRACCOEFS *frcPtrO; 
double longlist [ ] ; 

int nCross,nWindows,draw; 

/* Called by boxObserved. Gives observation times of box until t 

Assumes entry with satellite at equator, t at current time, 
longSat current longitude, and necessary parameters 

calculated 

in getConstants . */ 

^ int ilam, index; 

double longRel,tObs,frac,dlam; 

FRACCOEFS *frcPtr; 

DRAWDATA saw; 


SHOW "\n\nDay 


Hour 


Fraction") ; 
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while 

{ 


{ t < tLim ) 


longRel 


my fmod ( longSat - longLeft, TwoPi) ; 

/* Check sat. position relative to 


window */ . ^ T ^ rr 

if( longRel < windowLong && 

( (ilam - 

whichlntrvl (longRel, longlist, nCross, nWindows) 

* dlam = longRel - longlist [ilam] ; 

index - (nWindows — 1) ? 0 : ilam/ 8; 

tObs - t + tLeft [index] + dTdL[index] * (longRel 

- longlist [8*index] ) ; „ . 

frcPtr = frcPtrO + ilam; ^ * 

frac - frcPtr->bf + dlam * (frcPtr->mf + dlam 


frcPtr->qf) ; 
/*if ( draw ) 
getDrawData 


tellTF ( &tObs,4frac ); 

sdlam, drawPtrO+ilam, &saw) ; 

showPartSeen ( &saw ) ; 


} */ 

/* Go to next crossing of equator */ 
longSat - Normalong (longSat + crossLong) ; 
t +“ crossTime; 

} 

return ; 


void tellTF ( ptSec,pfrac ) 

/* Break tSec up into day, hour and send to output file. Later 
more 

processing may be done. */ 
double *ptSec, *pfrac; 

{ 

int day; 

double hour, intptr, fraction; 


day - ^ptSec / 86400. - (*ptSec < 
/* hour - (*ptSec - (double) day * 
hour - (*ptSec ) / 3600. ; 
fraction = modf(hour,&intptr); 
if (fraction > .5) 

hour - ceil (hour) ; 


0 .) ; 
86400.) 


/ 3600. 


else 

hour - floor (hour); 

/*if(hour " 24.) 
hour * 0 " ^ / 

SHOW "\n%2d’ %6.2f %6.3f",day,hour, *pfrac) ; 

/•SHOW "\t%6.3f", *pfrac) ;*/ 
return; 


} 


★ * 


★ * 

★ ★ 


* * 

★ * 


* ★ 
*/ 


i€ 


int NumObservO 



/* 

by 

{ 


Runs through specified box center latitudes and prints out the 
average number of times in 24 hours points in the box are seen 

satellite. */ 


extern int printLevel; 
extern double getFracO; 


double boxCtr,delLat,numObs, obsPerDay, latStepr latStep2; 

/* Print out some header information */ 
boxCtr - 0.; /* In degrees of latitude */ 

boxMatters - NUMOBS; /* Let setConstants we are running 
NumObserv */ 

printLevel « 1; „ , 

printf ("Working on latitude %7.2f\n ,boxCtr) ; 


if ( setConstants 0 ) 

return ( 1 ) ; 

delLat - phiN; /* Assumes box center is set at Equator in 

setConstants */ 

ObsPerDay = 24. * 3600./ tauN; 

SHOW "*\n") ; /* To help Cricket Graph ingestion */ 

SHOW "Latitude\tNumber\n") ; 
numObs - obsPerDay * getFracO; 

SHOW "%7.2f\t%7.3f\n",boxCtr, numObs) ; 

printLevel “ 0; 

latStep - 2.5 * 0.9999999; 

latStep2 - 1.0 • latStep; , ..o.. v 

for{ boxCtr ~ 2.5; boxCtr <“ 40.; boxCtr +~ latStep) 

^ if (boxCtr >23.) latStep = latStep2; 

printf ("Working on latitude %7 .2f\n", boxCtr) ; 

phiN = boxCtr*DtoR + delLat; 
phis - boxCtr *DtoR - delLat; 
if ( i “ boxCheckO ) 

{ 

if ( i !- 2 ) 

printf ("\nUnable to proceed because of 

box size."); 

break; 

numObs ” obsPerDay * getFracO; 

SHOW "% 7 ,2f\t%7 .3f\n", boxCtr, numObs) ; 

) 

return (0) ; 


* * * * * */ 


WhereSat ( ) 

t* Prints out latitudes and longitudes of satellite nadir, right 

swath edge, and time, starting satellite at equator and time t 
longitude 
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longSat, correcting for shift in longitude to put boxCtr at 0! 
long. 

{ 

extern int direction; 
extern double 

t , dt ^ tSt art , RtoD, alpha , bet a ^ sinlnc^ cosine, longSat , longBoxR; 
double phiEdge, longEdge, side, longSatShift ; 
double 

tShift,tSat,pSat,lSat,fabs {) ,sin() ,cos () ,atan2 () ,Normalong () ; 
double s, c, temp, tLimit ; 

boxMatters * WHERESAT; /* Tell setConstants that WhereSat is 
being run. */ 

setConstants () ; 

longSatShift * longSat; 
tShift » t; 


tSat * -travelTime; 

/* 

Gets the satellite back 

from equator to */ 


its starting position. 

ISat * 0 . ; 

/* 

"travelTime" is*/ 



pSat » 0 . ; 

!* 

left over from call to 

getlLongO in setC. */ 


Move to starting time. */ 

tSat +* tStart - t; 

/* 

SHOW 




"\nMins. \tlatNadir\tlongNadir\tlatR\tlongR\tlatL\tlongL\n") ; 

tLimit * tSat + tLim; /* Set tLimit so that satellite 


orbits that long*/ 

while ( tSat <- tLimit ) 

{ 

s = sin( (temp * tSat * alpha) ); 

c - cos ( temp ) / 

pSat * asin( sininc * s ); 

ISat * atan2 ( cosInc*s, c) ; 
if( fabs(lSat) > Pi2 ) 
direction * 0; 

else 

direction » 1; 

ISat “ Normalong( ISat - beta*tSat + longSatShift ); 
SHOW 

”%10 . If \t%7 .2f \t%7 . 2f \t" , (tSat+tShift) /60 . , pSat*RtoD, lSat*RtoD) ; 
side - 1 . ; 

SwathEdge{ pSat, ISat, side, &phiEdge, &longEdge ); 
longEdge - Normalong (longEdge) ; 

SHOW ”%7 .2f\t%7 .2f\t",phiEdge*RtoD, longEdge*RtoD) ; 
side - -1.; 

SwathEdge { pSat, ISat, side, &phiEdge, &longEdge ) ; 
longEdge = Normalong (longEdge) ; 

SHOW "%7 .2f\t%7 .2f\n",phiEdge*RtoD, longEdge*RtoD) ; 
tSat +* dt; 

} 

return; 

} 

* ★ * * * ★ ★ ★/ 


int BoxObservedO 

/* Runs Orbit Sat when called from main, 
values 

are specified in setConstants. */ 


Assumes all parameter 
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{ 

int nCross, nWindows t index^ gap/ vNum^ sLab^ dirx, bs, draw; 
double boxCtr, longlist [16] , temp, l,phiSw, longSw,pos; 

CROSSINGS *tableCross () , ♦pCross, *pCrossO; 

DRAWTEMPLATE ’^drawTempPtrO, *drwPtr; 

FRACCOEFS rcCoefPtrO, *frcPtr; 

boxMatters - BOXOBS; /* Inform setConstants that 
BoxObserved will be run */ 
printLevel - 1; 
draw - NO ; 

if ( setConstants 0 > 0 ) 
goto error; 

if ( (pCrossO » tableCross (inCross, AnWindows) ) “ NULL ) 
goto error; 

SHOW "\nTable of Equatorial Crossings. nCross=»%2d”, nCross) ; 
SHOW "\n\nindex longit. vertex side dirx”) ; 

for ( index*»0,pCross « pCrossO; index < nCross ; 
pCross++, index++ ) 

{ 

1 - pCross->longEq * RtoD; 
vNum “ pCross->vertex; 
sLab * pCross->SwEd; 
dirx « pCross->dir; 

SHOW "\n%3d %9.3f %4d %4d %4d", index, 1, vNum, sLab, dirx) ; 

} 

/* Allocate memory for tables 
drawTempPtrO -= (DRAWTEMPLATE *) malloc ( (unsigned) 

( (nCross-1) * 

sizeof 

( DRAWTEMPLATE ) ) ) ; 

if ( drawTempPtrO — NULL ) { 

printf ("\nTrouble getting memory allocation for 
drawTeitplate ! ”) ; 

goto error; 

) 

frcCoefPtrO - (FRACCOEFS *) malloc ( (unsigned) ( (nCross-1) 

★ 

sizeof 

( FRACCOEFS ) ) ) ; 

if ( frcCoefPtrO — NULL ) { 

printf ("\nTrouble getting memory allocation for 
fracCoefs ! ") ; 

goto error; 

) 

/* Construct tables of coefficients for running 

orbitSat */ 

if ( getCoefs ( 

nCross,nWindows,pCrossO,drawTempPtrO, frcCoefPtrO,draw ) ) 
goto error; 

/* Set up list of relative longitudes */ 
pCross « pCrossO; 
temp - pCrossO“>longEq; 

for ( index*0 ; index<nCross ; index++, pCross++ ) 

longlist [index] » my_fmod ( pCross->longEq - temp , 


TwoPi ) ; 


/* Get window size. Satellite w/ longitudes beyond 

this sees nothing. */ , , ^ . 

windowLong * longlist [ nCross-1 ] / 

/* Get time interpolation */ 
gap = (nWindows“=2) ? 7 : nCross-1; 

for°^^index=0 ; index<nWindows ; index++,pCross++ ) 

^ getLongO ( 0 . ,pCross->longEq, l-pCross->dir,pCross- 

>latTr)t, &temp ); 

tLeft [index] = travelTime; 

JetSongo" (^0^;pCross->longEq, l-pCross->dir,pCross- 
>latTrk., stemp^^ , ^ ^ „ (travelTime - tLeft [index] ) / 

longlist [gap] ; 

longLeft - pCrossO->longEq; 

if { free ( (void *) pCrossO ) ) /* Release memory used by 
pCross */ 

goto error; 

drawTempPtrO?frcCoefPtrO, longlist, nCross,nWindows, draw ) ; 

return (0) ; 
error: 

return (1) ; 


) 



Technical Report 


Section 2 


The feasibility of obtaining the longwave radiation budget 


the surface from satellite data. 



Reprinted from Journal OF Climate, Vol. 3, No. 12, December 1990 

American MeteorotopcaJ Soaery 


Relationship between the Longwave Cloud Radiative Forcing 
at the Surface and the Top of the Atmosphere 

Harshvardhan 

Department of Earth and Atmospheric Sciences, Purdue University. fVest Lafayette, Indiana 

David A. Randall and Donald a. Dazlich 

Department of Atmospheric Science, Colorado State University, Fort Collins. Colorado 
(Manuscript received 2 April 1990, in final form 20 June 1990) 

ABSTRACT 

Attempts to map the global longwave surface radiation budget from space have been thwarted by the presence 
of clouds. Unlike the shortwave, there is no physical relationship between the outgoing longwave and the surface 
longwave under cloudy skies. Therefore, there is no correlation between spatial and temporal averages of the 
outgoing longwave radiation and net longwave radiation at the surface. However, in regions where a panicular 
cloud regime exists preferentially, a relationship between the mean longwave cloud radiative forcing (CRF) at 
the top of the atmosphere and at the surface can be shown to exist. Results from a general circulation model 
suggest that this relationship for monthly means is coherent over fairly large geographical areas. For example, 
in tropical convective areas, the longwave CRF at the top is very large, but at the surface it is Quite small because 
of the high opacity of the lowest layers of the atmosphere. On the other hand, in areas of stratus over cool ocean 
surfaces, the longwave CRF at the top is very small but at the surface, it is quite substantial. 

To the extent that the cloudiness simulated in the model mimics the real atmosphere, it may be possible to 
estimate the monthly mean longwave CRF at the surface from the Earth Radiation Budget Experiment cloud 
forcing at the top. The net longwave radiation at the surface can then be mapped if monthly means of the dear- 
sky fluxes are obtained by some independent technique. 


1. Introduction 

Currently, there is a need in the atmospheric and 
oceanographic communities for reliable estimates of 
the components of the surface radiation budget ( WCP- 
92 1984). This is particularly crucial in the study of 
the coupled atmosphere-ocean system in the tropics 
(NAS 1983). Since there are very few surface stations 
measuring the radiation budget routinely and reliably, 
it has been realized that space-based observations are 
the only means to obtain global coverage. Several in- 
vestigators have had considerable success in obtaining 
the solar flux components at the surface using radiances 
measured at the top of the atmosphere by operational 
satellites (e.g., Gautier et al. 1980; Raschke and Preuss 
1979; Tarpley 1979; Pinker and Ewing 1985; Justus el 
al. 1986). This success can be explained on theoretical 
grounds. The atmosphere is a conservatively scattering 
medium over a major portion of the solar spectrum, 
and the total atmospheric absorption in the presence 
of clouds is not very different from the clear-sky value 
(Ramanathan 1986). This is because clouds absorb at 
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roughly the same near infrared wavelengths as water 
vapor, and this absorption is therefore at the expense 
of vapor absorption below the cloud deck. 

In the case of longwave radiation, it is not at all 
evident that the surface radiation budget components 
can be obtained from top of the atmosphere (TOA) 
radiance measurements. Even for clear-sky conditions, 
it may be shown that the bulk of the downward long- 
wave radiation emanates from the lowest hundreds of 
meters of the atmosphere (Schmetz et al. 1986; Gupta 
1989). Fairly accurate estimates of the near surface air 
temperature and humidity, and the sea-surface or 
ground temperature are prerequisites for the compu- 
lation of the longwave fluxes at the surface. This prob- 
lem is funher compounded by the presence of clouds. 
The downward emission from the cloud base is a major 
component of the downward longwave flux unless the 
cloud base is very high and the boundary layer is very 
moist. Since there is as yet no proven technique for 
locating cloud bases from space, investigators have been 
forced to rely on other means to estimate the downward 
emission. 

There is now considerable literature available on es- 
timates of the surface longwave fluxes, both regionally 
and globally. All of them rely on estimates of near sur- 
face conditions or the vertical profile of temperature 
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The theoretical basis for such a relationship in the 
case of solar radiation has been discussed by several 
authors, e.g., Gautier et al. ( 1980). All techniques for 
inferring the solar radiation at the surface rely on es- 
sentially the same principle. The chief uncertainty is 
in assigning absorption by water vapor in the clear at- 
mosphere and by clouds and vapor for the cloudy case. 
As pointed out by Ramanaihan ( 1986), compensating 
effects result in this absorption being fairly insensitive 
to location and in fact, cloudiness conditions. Modeling 
studies quoted by Ramanathan ( 1986) and Weare 
( 1989) confirm that there should be a strong correla- 
tion, at least for monthly means, between the net solar 
radiation at the top of the atmosphere, which is the 
absorption by the surface and atmosphere, and the net 
at the surface which is simply the surface absorption. 
However, the same modeling studies show that there 
is no correlation whatsoever between the net longwave 
radiation at the top of the atmosphere and the net 
longwave at the surface. 

Although Weare (1989) confirmed Pinker et ai.’s 
( 1985 ) finding that there is a strong correlation between 
the net total radiation at the top of the atmosphere and 
that at the surface, this is explained by the tact that the 
variation in the net total radiation at the surface is 
dominated by solar radiation. 

An illustrative way to show the correlation or lack 
thereof is with a scatter diagram of the radiation field 
at the top of the atmosphere plotted against the surface 
quantity. This can be done for contemporaneous mea- 
surements made at a surface site and satellite mea- 
surements of the same general area. This is, in fact, the 
usual manner in which regression equations for the 
surface radiative fluxes are calibrated. The regressions 
could also be for temporal means, or if several surface 
sites are available, an area average may be considered. 
Model simulations of the radiation fields can also be 
shown on a scatter plot. 

Following Ramanathan (1986), we show in Figs. 
la,b, the monthly mean net longwave radiation at the 
top and the surface as simulated by the UCLA/CSU 
General Circulation Model (GCM). Details of the ra- 
diation and cloud parameterization in the model and 
simulated fields of cloudiness and radiative fluxes may 
be found in Harshvardhan et al. ( 1989). A brief de- 
scription of the GCM may be found in the Appendix 
of Randall et al. ( 1989). Each point represents the 
mean for one grid point of the model for January and 
July, respectively. The horizontal resolution of the 
model is 4® in latitude by 5® in longitude. It is evident 
that the longwave radiation at the top is uncorrelated 
with that at the surface, and Ramanathan (1986) has 
shown that this holds true even if selected regions are 
considered in isolation. Although total fluxes are not 
correlated, the hope that window radiances provide in- 
formation on near-surface conditions has led to the 
development of hybrid techniques in which satellite 
measurements are used to reconstruct the temperature 


and humidity to compute clear sky fluxes using bulk 
formulae or a radiative transfer model and modifica- 
tions for observed or climatological cloudiness. Fung 
et al. ( 1984) reviewed several bulk formulae used for 
this purpose and compared them with a radiative 
transfer model. They concluded that uncertainties in 
the longwave fluxes due to lack of knowledge about 
cloud properties is larger than that due to expected 
uncenainties in temperature or humidity. Frouin et al. 
( 1988) compared several methods arranged in order 
of decreasing sophistication to compute the downward 
longwave surface fluxes for comparison with obser- 
vation taken during the Mixed Layer Dynamics Ex- 
periment (MILDEX). A technique that can be used 
only during daytime involving an estimate of the liquid 
water column amount in the cloud had the best cor- 
relation with observed fluxes. A relationship between 
cloud properties in the solar and cloud thickness was 
also used by Schmetz et al. ( 1986). Chou ( 1985) used 
cloud thicknesses based on London’s (1957) clima- 
tology whereas Darnell et al. ( 1986) and Gupta ( 1988 ) 
used a fixed cloud thickness of 50 mb. 

Apan from the problem of locating the cloud base, 
there is also the question of the fractional coverage of 
clouds that should be used in weighting the clear and 
cloudy sky contributions. Recently, Wu and Cheng 
( 1989) have used HIRS2/MSU sounding data to ex- 
tract all the relevant cloud parameters needed for their 
algorithm to compute the downward longwave flux 
globally for January and July 1979. There is no satis- 
factory verification for any of these methods. 

In this paper we present a theoretical framework that 
avoids the explicit computation of cloud fraction and 
the location of cloud base. Our hope is that global maps 
of the monthly mean net longwave flux at the surface 
may eventually be obtained in this manner. 

2. Surface radiation budget 

In general, the surface radiation budget has four 
components: the upward and downward directed 
longwave radiation, the incoming solar radiation at the 
surface, and the reflected solar radiation. The latter 
two components are, of course, present only during 
daylight hours. For some applications it is sufficient to 
know the net radiation, which is the algebraic sum of 
the upward and downward components. 

The net solar radiation at the surface is the surface 
absorption and varies from 0 to 1200 W m“^, while 
the net longwave at the surface is a radiative loss by 
the surface that can range from —200 W for a 
warm desert surface under clear-sky conditions to near 
zero in the presence of low thick clouds. The theory 
underlying techniques for estimating these quantities 
from space based observation systems is the relation- 
ship between the up welling radiation at the top of the 
atmosphere and the components of radiation at the 
surface. 



December 1990 


harsh VARDHAN, RANDALL AND DA2LICH 


1437 



-20 0 20 40 60 80 100 120 140 160 180 


Net LW at Surface (Wm’^) 



-20 0 20 40 60 80 100 120 140 160 180 


Net LW at Surface (Wm‘^) 

Fig . 1 . Scatter plots of the monthly mean net longwave radiation 
at the top of the atmosphere vs the net longwave radiation at the 
surface for (a) January, and (b) July. Results are from a simulation 
of the UCLA/CSU GCM and each data point on the plot represents 
a 4° Latitude X 5“ Longitude grid point. 


and humidity profile and a radiative transfer code is 
used to compute the downward longwave flux (e.g., 
Frouin et ai, 1988; Gupta 1989: Wu and Cheng 1989). 
As stated earlier, all techniques suffer from uncenain- 


ties or even total lack of knowledge concerning cloud 
cover and cloud base location. 

3. Cloud Radiative Forcing 

The role of clouds in modifying the radiation field 
can be studied through the analysis of the cloud radia- 
tive forcing (CRF). Although originally introduced in 
connection with the exiting shortwave and longwave 
fluxes at the top of the atmosphere (Charlock and Ra- 
manathan 1985; Ramanathan 1987), the cloud radia- 
tive forcing at the surface maybe defined as 

surface CRF = surface flux “ clear-sky surface flux. 

(i) 

In the above, the flux at the surface may be shortwave, 
longwave, or total, and could refer to the downward 
or net flux at the surface. The clear sky flux in a model 
may be obtained in one of two ways. The radiation 
fluxes may be computed at every grid point for the 
diagnosed cloud cover and also for zero cloud cover at 
the same time. The clear sky values may then be sub- 
tracted from the actual values to yield the CRF. This 
has been called ‘‘Method 11” by Cess and Potter ( 1987) 
and is the method used in this study. An alternate ap- 
proach is to sample the clear-sky fluxes only when 
clouds do not occur and build up a clear-sky clima- 
tology for a panicular grid point over an integrating 
time interval, such as a month. This is “Method I” 
and is used in the Earth Radiation Budget Experiment 
(ERBE; Ramanathan 1987) and is, of course, the only 
option in any measurement program. 

The importance of cloud radiative forcing at the 
surface may be illustrated by examining the energy 
budget at the atmosphere-ocean interface and impli- 
cations regarding the required meridional energy 
transport in the ocean for ^obal balance over the an- 
nual cycle. Esbensen and Kushnir ( 1981 ) have com- 
puted the individual components of the heat budget 
over the oceans using standard bulk formulae ( Sellers 
1965; Budyko 1974). Over a complete annual cycle, 
assuming that there is no heat storage in the oceans, 
an implied poleward energy transport across latitude 
belts maybe computed based on the excess or deficit 
of energy at each latitude. The results of Esbensen and 
Kushnir are shown in Fig. 2 by the solid line labeled 
EK. An independent estimate of the meridional oceanic 
transport has been made by Carissimo et al. ( 1985) 
who have updated the earlier work of Oort and Vonder 
Haar ( 1976). This is shown by the solid line in Fig. 2 
labeled COV. They computed the poleward transport 
of energy by the oceans and atmosphere using mea- 
surements of net radiation at the top of the atmosphere. 
The procedure for estimating the atmospheric com- 
ponent is described in Oort and Peixoto ( 1983). The 
oceanic transport is then obtained as a residual. 
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FiC . 2. Estimates of the annual mendional energy transport in the 
oceans. Solid line is from Carissimo et al. { COV ) and Esbcnsen and 
Kushnir ( EK). Dashed line is from an atmospheric general circulation 
model. Dotted line is from the same model but without cloud radiative 
forcing (CRF), Nonhward transport is positive. 

Figure 2 also shows by dashed lines the implied me- 
ridional transport in the oceans from a five-year annual 
cycle run of an atmospheric general circulation model 
with prescribed seasonally varying sea surface temper- 
atures. The procedure used to obtain the transport is 
the same as for the EK curve, except that the simulated 
components of the energy budget have been used. Al- 
though the magnitude of the transport in the Northern 
Hemisphere is within the range of the other estimates, 
it is the implied transport in the Southern Hemisphere 
that is truly striking. Over the annual cycle, the net 
oceanic transport is northward at all latitudes! The sit- 
uation is improved markedly by the removal of the 
radiative contribution of clouds, as shown by the dotted 
line. An examination of the radiation budget in the ^ 
Southern Hemisphere midlatitudes indicates that there J 
is a gross underestimation of cloud cover in that region. ^ 
Since the sea surface temperature is prescribed in the ^ 
model, these errors are of no consequence as far as the eS 
atmospheric simulation is concerned. However, errors ^ 
such as this could prove to be disastrous in a coupled ^ 
interactive atmosphere-ocean climate model, ^ 

Figure 3 is a scatter plot of the monthly mean short- 
wave CRF at the top of the atmosphere versus the 
shonwave CRF at the surface as simulated by the 
UCLA/CSU GCM. Each point represents the monthly 
mean (July in this case) value for each grid point. The 
CRF has been computed using Method II, in that dear- 
sky values were calculated and stored for every com- 
putation of the radiative fluxes. A linear correlation is 
evident and simply indicates that a loss of solar inso- 
lation at the surface due to cloud cover appears as an 


increase in the solar radiation reflected to space. The 
net shortwave at the top and surface vary in lock step 
with changing cloud cover. It is this feature that is pri- 
marily responsible for the success of surface shortwave 
insolation techniques based on satellite measurements 
of reflected radiances. 

It is instructive to investigate the simulated relation- 
ships between the longwave CRF at the top and surface, 
to see if any correlations are predicted by the model. 
The shortwave correlations could have been anticipated 
by inspection of Fig. 10 in Harshvardhan et al. ( 1989) 
in which the zonal means of the cloud radiative forcing 
from the same model are displayed. The zonal mean 
shortwave cloud radiative forcing at the top of the at- 
mosphere is virtually identical to the forcing at the sur- 
face with the atmospheric component alone being es- 
sentially zero. As explained in section Tthis is because 
the cloud layer absorbs solar radiation at the expense 
of water-vapor absorption below the cloud deck and 
the vertically integrated column absorption is therefore 
virtually unchanged in the presence of clouds ( although 
the heating profile is certainly quite different). 

Figure 4 shows a scatter plot of the monthly means 
of the longwave CRF at the top of the atmosphere ver- 
sus the surface CRF for every grid point in the GCM. 
Note the radically different distribution from Fig. 1, 
which is a plot of the net longwave radiation. Although 
the scatter is incoherent as the shortwave CRF shown 
in Fig, 3, there are unmistakable correlations among 
groups of points. It appears, in fact, that certain clusters 
of points are linearly organized, albeit with different 
slopes. These groups of points represent grid points at 
which, for the month considered, a preferred type of 
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Fig. 3. Scatter plot of the July mean shonwave cloud radiative 
forcing (SVVCRF) at the top of the atmosphere vs the SWCRF at 
the surtace as simulated by the UCL^/CSU GCM for all gnd points. 
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Figure 5a shows the scatter plot for two zonal st^s 
in the tropics in July. The points are all the GCM 
centered at 10° and 14°N (144 in 
the same plot for two bands centered at 58 and bl b. 
The former includes areas of high convective clouds 
whereas the latter is a region of low stratus. The long- 
wave CRF at the top of the atmosphere in the tropics 
can be very large since the cold cloud tops diminish 
the outgoing longwave radiation considerably. How- 
ever, at the surface, in the tropics, clouds have a min- 
imal effect on longwave radiation because the lowest 
layers of the atmosphere have a very high water-vapor 
content, essentiaUy saturating most of the longwave 
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Fig. 4. As in Fig. 3 but for the longwave cloud radiative forcing 
(LWCRF) for (a) January and (b) July. 


cloud cover exists such as low or convective. This be- 
comes dear when subsets ot the global grid are consid- 
ered where one would expect this sort of preference to 
hold. 
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Fir 5 As in Fig 4 for Julv but only for gnd points in zonal 
centered at (a) 10° and 14°N and (b) 58° and 62°S. 
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spectrum for downwelling radiation. The presence of 
an emitting cloud base does not result in much en- 
hanced downward longwave radiation. Moreover, in 
the model, low-leve! convective clouds are ignored in 
that the cloud fraction is set to zero, so the cloud base 
is quite high unless stratiform clouds occur simulta- 
neously. The situation at high latitudes in drier con- 
ditions is just the opposite. Stratus clouds do not modify 
the outgoing longwave radiation much, but radically 
alter the downward longwave flux towards the surface. 
Figure 5b also indicates that middle and high clouds 
occurred infrequently, if at all, during this time period 
at the grid points considered. 

The points in both Fig. 5a, b are not clustered at the 
upper right-hand comer of the distribution, but are 
spread out fairly even along the line of correlation. The 
position of a point is an indication of the monthly mean 
cloud cover for the grid point in question. From Eq. 

( 1 ) we note that the first term on the rhs includes the 
cloud cover implicitly, such that if there is zero cloud 
cover at a grid point, the CRF at both the top of the 
atmosphere and surface will be zero. 

The character of the correlations shown in Fig. 4a,b 
can be best appreciated by considering the simulated 
longwave CRF at the top of the atmosphere and surface 
for standard atmospheric profiles in which clouds are 
inserted at different levels in isolation. We have chosen 
to show standard profiles from McClatchey et al. 

( 1972) and cloud layers corresponding to the 9-layer 
UCLA/CSU GCM. Optical properties of the clouds 
are as in the GCM ( Harshvardhan et al. 1989); essen- 



LW CRF at Surface (Wm*2) 

Fig. 6. The locus of points representing the variation of LWCRF 
at the top of the atmosphere vs the LWCRF at the surface when 
cloud cover varies from 0% to 100% in each layer of a model in 
isolation. A midlatitude summer atmospheric profile is used and the 
cloud properties and layer positions are as in the UCLA/CSU GCM. 



Fig. 7. As in Fig. 6 but for a subarctic winter atmospheric profile. 


tially, the longwave emittance is a function of temper- 
ature and cloud thickness, but, for all practical pur- 
poses, this dependence exists only for temperatures well 
below freezing. For water clouds and ice clouds that 
form as convective anvils, the emittance is unity. 

Figure 6 is a plot of the longwave CRF at the lop of 
the atmosphere versus the surface CRF. It shows the 
locus of all points representing cloud cover in one layer 
of the model for a midlatitude summer atmospheric 
profile. There are eight different lines, one for each 
layer in which clouds are permitted to occur in the 
GCM. The origin represents the case of no cloud cover 
while the terminus of each line represents 100% cloud 
cover. The line with the lowest values of longwave CRF 
at the top and highest values of longwave CRF at the 
surface corresponds to a cloud in the lowest layer of 
the model. Clouds in higher layers proceed in a coun- 
terclockwise sense up to the topmost cloud layer in the 
model. 

The CRF of a cloud in the highest layer, with a top 
at 100 mb, is barely perceptible because the emittance 
is close to zero and even complete cloud cover hardly 
affects the outgoing longwave radiation. It should be 
stressed that points on a scatter plot would follow these 
lines only if cloud existence occurred in that panicuiar 
layer alone, not if clouds existed simultaneously in 
more than one layer. 

Figure 7 is similar to Fig. 6, but for a subarctic winter 
atmospheric profile. Figure 8 is for the particular case 
of tropical convective anvils, which in the model extend 
from about 400 mb to the top of the detrainment layer 
in the convective parameterization (Randall et ai. 
1989). These anvils then have the same base level, but 
have varying tops up to 100 mb. Low-level convective 
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Fig. 8. As in Fig. 6 but for convective anvils 
and a tropical atmospheric profile. 


clouds are completely ignored, at present, in the ra- 
diation computations. 

The simulations shown in Figs. 6-8 explain the na- 
ture of the scatter plots shown in Figs. 4, 5. The model 
tends to generate a preferred type of cloud at a partic- 
ular grid point during a particular season. Although 
simultaneous occurrence of multilayer clouds is sim- 
ulated, there are several regions in which a particular 
cloud type predominates. For example. Fig. 5a clearly 
shows the effects of tropical anvils. Some grid points 
along these zonal bands exhibit a slightly different cloud 
partem, but there is a large area of coherence. The 
same is true of Fig. 5b, which is dominated by low- 
level clouds, and may be compared with the lowest line 
in Fig. 6. 

Regions with a preferred type of cloudiness were 
identified by Hartmann and Short (1980) using the 
day-to-day variability of the outgoing longwave radia- 
tion and shortwave albedo. They found that the char- 
acteristic slope of the line of regression of changes in 
outgoing longwave radiation and shortwave albedo was 
a very sensitive indicator of low cloud regimes. 

IVe have shown that the characteristic slope of the 
line of regression of the longwave CRF at the top of the 
atmosphere and the surface CRF also identifies cloud 
regimes. Following Hartmann and Short ( 1980) maps 
of the ratio of the longwave CRF at the top of the 
atmosphere to the longwave CRF at the surface for 
January and July are shown in Figs. 9a,b, respectively. 
The map was produced from the data in the scatter 
plots shown in Figs. 4a,b. Regions for which this ratio 
is small can be identified with areas of predominantly 
low-level clouds or clouds in high latitudes in the win- 


ter. Large values of the ratio indicate convective anvils. 
The maps in Figs, 9a,b are thus an indication of the 
regional cloud climatology of the model. This identi- 
fication can perhaps be utilized to estimate the monthly 
mean net longwave radiation at the surface. 

4. Conclusions 

All methods for extracting the surface longwave ra- 
diation suffer from uncertainties regarding cloud cover 
and cloud base temperature. The analysis presented in 
this study suggests a technique that may be able to 
avoid this problem. The longwave cloud radiative 
forcing ( LWCRF) at the top of the atmosphere is cur- 
rently being compiled by the Earth Radiation Budget 
Experiment (ERBE). The mean downward longwave 
radiation can be computed for clear sky conditions if 
a mean profile of temperature and water-vapor mixing 
ratio is available. The upward longwave radiation at 
the surface may be computed from a retrieved ground 
temperature. The LWCRF at the top contains infor- 
mation on the mean cloud cover and cloud type for 
the time period considered. If monthly means are used, 
the LWCRF at the surface may be obtained from a 
global map of the LWCRF at the top if maps of the 
ratio, such as shown in Fig. 9, were available. This 
would not require explicit knowledge of the mean cloud 
fraction. Of course, acceptance of the results of a general 
circulation model is not a viable option. To some extent 
the standard deviation of daily mean values of the ratio 
of cloud forcing from the model could be used to judge 
the appropriateness of this technique for a particular 
grid point. This will involve reliance on the model that 
is probably not warranted at this stage. Some infor- 
mation can also be obtained, however, from the results 
of recent or planned field experiments. 

Certain elements required for this procedure to be 
tested are already in place. The data products released 
by the International Satellite Cloud Climatology Proj- 
ect ( ISCCP, Schiffer and Rossow 1983 ) can be directly 
used in this scheme. ISCCP provides a daily profile of 
temperature and precipi table water as well as surface 
temperature at a 2.5° X 2.5° horizontal resolution that 
may be used to generate clear sky downward longwave 
fluxes and upward fluxes, at least over the oceans where 
the once-a-day sampling is acceptable. As part of 
ISCCP, there have been field experiments to study cir- 
rus and marine stratus cloud regions (Starr 1987; Al- 
brecht et al. 1988), .At least for these cases, an estimate 
of the mean LWCRF at the surface for an integrating 
lime period comparable to ERBE can be obtained. On 
a global scale, of course, model simulations are the 
only choice. 

Since current needs for surface radiation budget es- 
timates are being determined primarily for application 
in the tropical oceans (WCP-92 1984), a start could 
perhaps be made to generate monthly mean surface 
longwave radiation maps for the tropical Pacific using 
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Fig. 9. Global maps of the ratio of the LWCRF at the top and at the surface for (a) January, 
and (b) July constructed from a simulation of the UCLA/CSU GCM. The ratio was obtained 
from points plotted in Fig. 4. 


the procedure outlined in this study. The method 
should have a greater chance of success in that region 
based on the correlation shown in Fig. 5a. A test could 
be made as soon as ERBE and ISCCP data are available 
for the same time period. 


To summarize, although it appears at first glance 
that longwave fluxes measured at the top of the at- 
mosphere provide very little information on surface 
fluxes, it is possible to extract information from an 
analysis of the mean cloud radiative forcing that can 
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be used to generate regional, if not global, maps of the 
net surface longwave radiation. We have not attempted 
to estimate quantitatively, the accuracy of the final 
product, but have merely suggested a methodology that 
utilizes cloud and radiation datasets that are currently 
being compiled under the auspices of ISCCP and 
ERBE. 
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Summary 


Global maps of the monthly mean net upward longwave radiation flux at the ocean 
surface have been obtained for April, July, October 1985 and January 1986. These maps were 
produced by blending information obtained from a combination of general circulation model 
cloud radiative forcing fields, the top-of-the-atmosphere cloud radiative forcing from ERBE and 
TOYS profiles and sea surface temperature on ISCCP Cl tapes. The fields are compatible with 
known meteorological regimes of atmospheric water vapor content and cloudiness. There is a 
vast area of high net upward longwave radiation flux (> 80 W m'^) in the eastern Pacific Ocean 
throughout most of the year. Areas of low net upward longwave radiation flu.x (<40 W m -) are 
the tropical convective regions and extra tropical regions that tend to have persistent low cloud 
cover. The technique used in this study relies on GCM simulations and so is subject to some ot 
the uncertainties associated with the model. However, all input information regarding 
temperature, moisture and cloud cover is from satellite data having near global coverage. This 
feature of the procedure alone warrants its consideration for further use in compiling global maps 
of the net longwave radiation at the surface over the oceans. 



1 


1. Introduction 

The surface radiation budget, i.e., the net solar radiation absorbed minus the net 
longwave radiation emitted, and its spatial and temporal variations are key parameters in climate 
and weather studies. This budget plays a major role in determining radiative heating, as well as 
sensible and latent heat fluxes over ocean and land surfaces. As a result, the net radiative flux 
constitutes an important boundary forcing for the general ocean circulation and a crucial 
parameter for determining meridional oceanic heat transport, ocean-atmosphere interaction and 
land-atmosphere interaction. Moreover, it is a useful parameter when addressing issues related 
to climate change due to CO 2 and other trace gases, and in the validation of radiation schemes 
used in climate models. Therefore, it is understandable for the atmospheric and oceanic 
communities to need reliable estimates of the surface radiation budget (WCP-92, 1984). 

Direct high-quality radiation measurements at the surtace are difficult to make, 
particularly over the oceans which cover more than 60% of the Earth's surface. Actually there 
are very few surface stations measuring the radiation budget routinely and reliably because ot the 
requirement of careful instrument calibration and temperature correction for the radiation, 
especially longwave, measurement. In addition, because of operating costs, it is not feasible to 
maintain a network of surface stations over the oceans. Although an attempt has been made to 
use ships to observe some meteorological parameters such as sea surface temperature, air 
temperature, specific humidity near the surface and even fractional cloud coverage, there has not 
been much progress in the measurement of surface radiation. Few ships measure radiation 
quantities because of special needs that require a dedicated tacility for this purpose. Moreover, 
regular ship observations are limited along commercial shipping lanes, and vast geographic gaps 
still exist, especially in the southern hemisphere. Consequently, the empirical formulas used to 
derive budgets have been validated only over limited regions, and when applied globally, large 
errors are inevitable. Therefore, direct measurement of shortwave and longwave radiation fluxes 
at the surface globally has not been possible. 

Since there are difficulties in obtaining radiative data from surface stations routinely and 
reliably, it has been realized that space based observations are the only means to have global 
coverage. However, because ot the intervening atmosphere, the surtace radiation budget is 
difficult to measure from satellites, whereas the top-of-the-atmosphere (TOA) radiation balance 
can be measured directly. Over the past decades, considerable effort has been expended in the 
tilobal measurement of the TOA radiation budget. Attempts at inferring the surtace radiation 
budget from space based measurements have only begun recently. 

There has been some success in obtaining the solar radiation budget at the surtace (e.g. 
Raschke and Preuss, 1979; Tarpley, 1979; Gautier et al.. 1980; Pinker and Ewing, 1985; Justus et 
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al., 1986). The progress in surface longwave radiation budget measurements trom space, 
however, has been much slower. Currently, some techniques are available to estimate the 
downward longwave component of the surface flux (Darnell el al.. 1983. 1986; Chou. 1985; 
Schmetz et al., 1986; Frouin et al., 1988; Ardanuy et al., 1989; Wu and Cheng. 1989; Breon et 
al , 1991). The net longwave component at the surface can be estimated by the ditlerence 
between the upwelling and downwelling fluxes. The upward component is determined directly 
from sea surface temperature since the oceanic surface emits essentially as a blackbody. The 
downward flux, however, is more difficult to obtain since it depends on many meteorological 
parameters such as atmospheric moisture, temperature and cloud cover. Because of the 
uncertainty of the measurement of these meteorological parameters, at present there is a need for 
improvement in the estimation of net longwave radiation fluxes at the surtace. . 

Two types of methods have been used to estimate the downward longwave flux at the 
surface; statistical and physical. As the name implies, statistical methods rely on correlations 
between fluxes and observed meteorological parameters. The physical techniques are based on 
modeling radiative processes occurring in the atmosphere (clear and cloudy atmosphere). The 
downward flux is computed from radiative transfer models which utilize parameters obtained 
from satellite radiance data. These parameters include temperature and water vapor mixing ratio 
profiles, fraction of cloud coverage and cloud emittance. However, all physical methods 
currendy under consideration have to make certain assumptions regarding both the presence of 
clouds and their vertical extent. Recent examples of these attempts are Chou (1985), Schmetz et 
al. (1986), Darnell et al. (1986), Gupta (1989), and Wu and Cheng (1989). The treatment ot 
longwave radiation transmittance in the presence of clouds becomes more complex since 
knowledge of cloud top and base heights and emittances are required. For this reason, it is 
important to determine the vertical profile of cloudiness as well as the horizontal distribution ot 
clouds and associated emittances. Unfortunately, determination of the vertical profile of 
cloudiness from space based measurements is difficult since overlap of cloudy layers is common 
in the real atmosphere. Therefore, because of the uncertaindes in assumed cloudiness, all these 
methods often give unreliable results. 

The method used here to obtain monthly mean quandties avoids the explicit computation 
of cloud fraction and the location of cloud base in estimating the downward longwave radiation 
globally (Harshvardhan et al., 1990). An advantage of this technique is that no independent 
knowledge or assumptions regarding cloud cover for a particular month are required. The only 
informadon required is a relationship between the cloud radiative forcing (CRF) at the top ot the 
atmosphere and that at the surface, which is obtained trom a general circuladon model (GCM) 
simulation. 
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Cl data provides a daily profile of temperature and precipitable water as well as surface 
temperature at a 2.5° x 2.5° horizontal resolution. This information from TOYS (TIROS 
Operational Vertical Sounder) is used to generate clear sky downward longwave fluxes globally 
and upward fluxes only over the ocean where the once-a-day sampling is acceptable. The 
diurnal cycle ol surface temperature precludes using this technique over land. The radiation 
code used to compute the clear sky downward flux is the one used in the UCLA/GLA (now 
CSU) GCM (Harshvardhan et al., 1989; Randall et al., 1989, 1991). 

A flow diagram of the technique used to obtain surface longwave fluxes is shown in Fig. 
1. Calculations start from the relationship between the longwave CRF at the top of the 
atmosphere and at the surface (Harshvardhan et al., 1990). Results for the months of January, 
April, July and October are used in this study and represent the simulated monthly characteristics 
of this relation over the annual cycle. The longwave CRF at the top of the atmosphere as 
obtained from FREE is combined with the relationship to obtain the longwave CRF at the 
surface for each of the four months. Then, by means of the radiative transfer model, with input 
meteorological parameters such as temperature profile and water vapor mixing ratio profile and 
sea surface temperature from ISCCP data as well as standard ozone vertical distributions, the 
clear sky downward longwave radiative flux and the clear-sky net longwave upward radiation 
flux at the surface are obtained. Because ozone is primarily confined to the stratosphere, its 
contribution to downward longwave flux is much less than that of other parameters such as water 
vapor in the lower atmosphere. Therefore, use of a standard ozone profile is justified. Trace 
gases such as methane and freons are not included in the code. This omission leads to an 
overestimate of the outgoing longwave clear sky flux by 5 to 10 W m*^ (Brieeleb, 1992). 

However, the downward clear sky flux is dominated by water vapor emission and the effect of 
trace gases is neglected here. 

In this study, the clear sky fluxes are computed for the atmospheric structure under 
cloudy conditions, but assuming a cloud-free sky. This is referred to as "Method 0" by Cess and 
Potter (1987). The ERBE clear sky climatology is constructed by averaging fields that are 
classified as being clear. The method of determining the samples to be included in this category 
introduces errors. Hanmann and Doelling (1991) have shown that fluctuations in temperature 
and water vapor can result in misclassification. As a result of this, the clear sky OLR over the 
oceans could be an overestimate by -3 W m'2 (Harrison et al., 1990). This bias will also appear 
in the surface CRF computed using the procedure outlined in Fig. 1. Once the CRF at the 
surface and the clear-sky radiative fluxes at the surface are available, based on the definition of 

the CRF shown in eq. (1), the actual radiative fluxes can be calculated simply as the clear-sky 
flux plus the corresponding CRF. 
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3. Results 

3.1 Longwave CRF at the Surface (global) 

The mean monthly distributions of the longwave cloud radiative forcing (CRF) at the 
surface are presented in Fig. 2 for the months of April, July and October 1985 and January 1986. 
They are derived by combining the longwave CRF at the top of the atmosphere from ERBE 
(Harrison et al„ 1990) and the ratios of the CRF at the top to that at the surface. Based on the 
definition in eq. (1), the longwave CRF is the difference between the mean longwave flux and 
clear sky longwave flux for the same period. Here the longwave CRFs at the top ot the 
atmosphere are retrieved from ERBE satellite data according to this definition. In the original 
ERBE data, there are a few regions in the tropics where clear sky longwave radiative fluxes at 
the top of the atmosphere are unavailable because of the lack ot clear sky pixels during the 
experiment period. Therefore, an interpolation is performed to make up the missing data using 
neighboring points. 

There are several noteworthy features of these distributions. First, regions with small 
surface longwave CRF are concentrated in the tropical and the subtropical oceanic areas (Fig. 2). 
Surface CRF is small in the tropics because the boundary layer there is moist and radiatively 
opaque even for clear skies. In the central Pacific Ocean, the areas with surface longwave CRF 
below 20 W m‘2 dominate throughout the year. This is also the case for the nonhem and central 
Indian Ocean in both April 1985 and January 1986. Areas with large surface longwave CRF 
occur over oceanic areas southwest of Indonesia, with a maximum of more than 60 W m “ in 
October and more than 80 W m'2 in July 1985. Both of these maxima correspond to areas of 
tropical convective activity. In the central Atlantic Ocean, the surface longwave CRF is below 
30 W m-2. 

Larger values of surface CRFs are found over the midlatitude continents than over oceans 
during October, January and April. This may seem surprising since cloud cover over the oceans 
is greater than that over land (Rossow and Schiffer, 1991) and the longwave CRF at the top of 
the atmosphere is larger over the oceans (Harrison et al., 1990). The results could be an artifact 
of the method but are not inconsistent with meteorological conditions for the following reason. 

The surface LWCRF defined in eq. (1) can be rephrased as: 

Surface LWCRF = CfFs^^c - Fsfclr) 

where C is the cloud fraction, F^^vc' is the surface longwave Oux for overcast conditions and 
Fsfclr is the corresponding value for clear skies. The surface LWCRF is thus the product ot two 
terms. The land-ocean difference in the quantity is therefore: 
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A LWCRF = CoCFjfovc ~ ^sfclr)© “ (^sfovc Psfclr)b 
where the suffix o refers to ocean and 1 to land. 

There is no doubt that Cq > Cj but the difference in the downward surtace longwave flux 
between clear and overcast conditions is larger over dry, cold continents than over moist, warm 
oceans (Chou, 1989). In fact, the maximum LWCRF at the surface occurs over the Himalayas 
(Wu and Cheng, 1989). Whether the larger flux difference fully compensates for the difference 
in cloud cover cannot be known without observational verification. Comparison with other 
modeled values is not sufficient grounds for rejecting this possibility. 

An example of the possibility mentioned above is a comparison of the net flux difference 
between overcast and clear skies for two disparate areas in the month of January 1986. We have 
considered an oceanic area in the western Pacific, which is the site of the TOGA-COARE 
experiment (WCP-92, 1984) and a continental area in eastern North America. Net upward 
longwave fluxes at the surface have been calculated using ECMWF (European Centre tor 
Medium Range Forecasts) analyses and the radiation code used in this study for clear sky 
conditions and two separate overcast cases, one with the cloud base at 850 mb and the other with 

the base at 931 mb. 

The clear sky net upward flux at the surface in the TOGA-COARE area is about 60 W 
m-2, while in the continental region in North America it is about 100 W m'^. The mean 
difference between overcast and clear in the oceanic area is 41 W m'2 for the 850 mb cloud base 
and 48 W m'2 when the base is at 931 mb. The corresponding mean differences in the 
continental area are 96 W m’2 and 95 W m'2 respectively. It should be noted that there is a 
pronounced temperature inversion in the analyzed fields for this region. These values are the 
terms in parenthesis in eq. (3), such that even if Cq is substantially greater than Cj, the 
forcing over land could exceed that over the oceans. 

There are two continents, Africa and Australia, which always have small surface CRFs 
throughout the year because of their large desert and semi-arid areas. Here the cloud fraction is 
low and the first term in eq. (2) determines the forcing. Fig. 2 shows the surface CRFs there to 
be often below 30 W m'^, especially in northern Africa, where the values are always below 20 W 
m-2. It is worthwhile noting that a persistent large surface CRF region exists in the tropical 
region of South America. Obviously, this results from the persistent cloud cover in this region. 
Over high latitude regions, some areas with the largest surtace longwave CRFs are found over 
the region poleward of 65°S, with a maximum of 150 W m'2 in July and some areas over the 
Arctic, with a maximum of more than 140 W m'2 in October and January. Two explanations are 
possible for this phenomenon. Physically, it is reasonable that persistent low clouds over these 
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regions lead to high surface longwave CRF. On the other hand, an underestimate of high cloud 
over the polar regions in the GCM simulation can also result in the surface longwave CRFs 
being somewhat unreasonably high. Also, FREE cloud forcing for these high latitudes may not 
be reliable. Finally, areas with the average range of 30 - 50 W m'- dominate most mid-latitude 
oceanic regions of both hemispheres. This is attributable to the increase in oceanic stratus in 
these regions. 

3.2 Downward Longwave Radiation Flux at the Surface (global) 

Downward longwave radiative fluxes at the surface for cloud-free skies, as calculated by 
the radiation code in the UCLA/GLA GCM in conjunction with the input meteorological 
parameters from the ISCCP satellite data and U. S. standard atmosphere (COESA, 1976), are 
presented in Fig. 3 for April, July and October 1985 and January 1986 respectively. 

Over the oceanic areas, clear sky surface downward longwave fluxes have pronounced 
zonal distributions, especially in mid-latitudes and near polar regions. In the tropical and 
subtropical areas, the regions with large downward fluxes (larger than 400 W m'-) are centered 
over Southeast Asia in April and October 1985 and shift a little northward in July 1985 and, as 
expected, a little southward in January 1986. This is not surprising since the surface downward 
longwave fluxes for clear skies are related closely with seasonal changes in temperature and 
moisture. In the northern spring (April) and fall (October), for the quite symmetrical solar 
irradiance distribution at that time, surface downward longwave fluxes also have a symmetrical 
distribution with respect to the equator. In the northern summer (July) and winter (January), the 
maps of clear sky surface downward fluxes show a shift northward and southward respectively 
with the changing seasons. This result is consistent with the fact that clear sky downward fluxes 
are determined by the near surface temperature which is closely related to the incident solar 
energy. Evidently, the high values (larger than 400 W m‘^) in Southeast Asia correspond to 
areas where, due to the high surface temperature of islands, the atmosphere is warmer than 
elsewhere along the same latitudinal belt. In particular, there >s an area with values larger than 
420 W m*2 centered east of Papua New Guinea in January 1986, corresponding to very warm 
and moist near surface conditions. In contrast, the low values (less than 200 W m"^) near the 
polar regions correspond to areas where the atmosphere is cold and dry throughout the year. 

Over the continents, symmetrical distributions of clear sky downward fluxes disappear 
and more complicated features of distributions emerge, corresponding to the land surface. There 
is a pronounced center of low values of downward fluxes over the Tibetan Plateau in Asia 
throughout the year because the mean elevation of the Tibetan Plateau is more than 4000 m 
above sea level so that surface air temperature is quite low except in summer. However even 
though high temperatures often occur there in summer, its high elevation makes the air still very 
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dry. Since the major contribution to longwave downward Hux at the surface comes from the 
water vapor in the lower atmosphere, the dry air over that region results in an area of low values 
of downward fluxes. Over southern Africa, relatively higher values along the same latitude 
result from the higher air temperature and more moist air. In nortliem Africa, although the 
surface temperature is very high, the dry air, due primarily to its desert areas, results in 
downward fluxes that are lower along the same latitude. In addition, a pronounced trough line 
along the Rocky Mountains in Fig. 3 represents higher altitudes and a drier atmosphere in that 
region compared to surrounding areas. 

Surface downward longwave fluxes, as derived from the formula described in eq. (1) 
with given clear sky fluxes and the longwave CRF at the surface, are presented in Fig. 4 for 
April, July and October 1985 and January 1986 respectively. Over the oceanic areas, the zonal 
distribution for the clear sky case is no longer apparent. The presence of clouds increases the 
downward fluxes at the surface globally throughout the year. Fluxes with values larger than 400 
W m"“ are found over areas associated with the intense convective cloud systems, such as the 
ITCZ, as well as the summer and winter monsoon areas. It is reasonable because, in these 
regions, clouds occur frequently and the air temperatures are high. An area with values larger 
than 440 W m‘- is found over northeast Australia in January 1986, corresponding to the high 
temperatures in the austral summer. Also, an area with high values is centered over southwest 
Indonesia in July 1985. 

Compared with clear sky maps, larger changes occur over land than over the oceans 
throughout the year. In addition, the regions of high values are still distinguishable in these 
months except in April 1985. Low downward fluxes (less than 280 W m'-) are found over high 
latitudes and polar regions, where the precipitable water is low and air temperature is cold. As 
seen in Fig. 4, contour lines of downward fluxes are spaced very densely in high latitudes during 
the whole year. As we mentioned before, these sharp changes perhaps result from uncertainties 
in modeling the meteorological parameters in these regions, an unavoidable consequence of the 
inability to distinguish cloud cover from background snow and ice and hence compute the cloud 
forcing. In addition, seasonal variations of surface downward fluxes are still notable in Fig. 4, 
even though they are not as remarkable as in the clear sky case. 

3.3 Net Upward Longwave Radiation Flux at the Surface (ocean only) 

The clear sky net upward longwave flux at the surface is the difference of the upward 
flux minus the downward surface flux assuming no clouds. The upward flux here is computed 
from the sea surface temperature, which is a reponed parameter in ISCCP Cl satellite data. 
There are three surface temperatures (TS) from imaging radiometers provided by ISCCP: mean 
TS from a clear sky composite, mean TS for IR-clear pixels, and mean TS for VIS/IR-clear 
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pixels for only day time. Also listed is a TOYS surface temperature. In order to compute the 
surface emission, it is not necessary to use any ot these fields but instead rely on an independent 
source such as the sea surface temperature provided by NOAA’s Climate Analysis Center 
(CAC). This is a blend of in situ data. Advanced Very High Resolution Radiometer (AVHRR) 
satellite data, and ice data (Reynolds, 1988). Fig. 5 shows the zonally averaged sea surface 
temperatures for January 1986 for the three ISCCP fields and the CAC field. In high latitude 
regions, the CAC sea surface temperature is much higher than the ISCCP values because ice 
surface temperature is set to be the freezing point of sea water (-1.8°C). Fig. 5 indicates that the 
clear sky composite TS is warmest among the three ISCCP sea surface temperatures, while the 
IR clear sky TS is the coldest. This is to be expected since IR clear pixels are contaminated by 
low level clouds. We chose the VIS/IR clear sky as the most representative quantity for this 
study. 

The difference between the surface emission computed from daydme mean TS for 
VIS/IR clear pixels and the CAC sea surface temperature is presented in Fig. 6 for January 1986. 
This difference is a measure of the uncertainty one may expect in the upward longwave radiation 
flux computed using different sources for the sea surface temperatures. The dominant feature of 
this map is a set of biases of less than 10 W m'^ present over the tropical and subtropical regions 
except around southeast Asia, where the surface emission from mean TS for VIS/IR clear pixels 
exceeds that from CAC temperature by as much as 30 W m'^. This could be a result of 
differences in the particular algorithms used by ISCCP and CAC to account for water vapor 
absorption in the atmospheric window and cloud screening procedures. Also any ship 
observations used are not necessarily representative of the skin temperature of the ocean surface 
but this effect should be small. The areas with differences larger than 10 W m'^ are found over 
high latitudes in the Pacific and Atlantic Ocean. Positive differences (less than 10 W m'^) are 
located over the central Pacific ocean and eastern Indian ocean, whereas areas with negative 
difference are found in other locations. The maximum absolute difference is around 5%. 
However, it is disturbing that the uncertainty in the upward component of the surface longwave 
flux is of the order ot the accuracy in the net radiation requested by the scientific community 
(WCP-92, 1984). 

Fig. 7 shows the clear sky net upward longwave fluxes at the surface for April, July and 
October 1985 and January 1986 respectively. Physically the maps of clear sky net upward 
longwave fluxes primarily reflect the distribution of water vapor content in the boundary layer. 
The area with lowest values is found over Southeastern Asia throughout the year. In this region, 
even in the absence ot clouds, the clear sky net upward longwave fluxes are quite small because 
of the high water vapor mixing ratio near the surface. Areas with values larger than 100 W m‘2 
are found over the mid-latitude subsidence zones between 15° and 45° latitude in both 
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hemispheres in April, October 1985 and January 1986, and in the southern hemisphere in July 
1985. These are regions in which a dry atmosphere overlies a moderately warm ocean surface. 
Also, in July 1985, there is a minimum of less than 40 W m'^ off the coast of central America 
and, in the other three months, a minimum ot less than 40 W m‘~ just along the west coast of 
Colombia, Ecuador and Peru; all these minima correspond to the presence of a relatively moist 
atmosphere. 

Net upward longwave flux, a difference field of the clear sky net upward flux minus the 
surface longwave CRF (Fig. 2) as described in eq. (1), is presented in Fig. 8 for April, July and 
October 1985 and January 1986. In this case, the areas with smaller values are still found over 
Southeast Asia. Here, the near saturation of the boundary layer results in low values of the net 
upward longwave flux for both clear and cloudy conditions resulting in low values in the 
monthly mean. Also, small net upward fluxes with values less than 40 W m'^ are found in high 
latitudes and polar regions where there is persistent cloud cover and temperatures are low 
throughout the year. Values larger than 80 W m'^ are found over some areas in the tropical and 
mid-latitude zone. It is worth noting that the negative values near the Antarctic in July 1985 
result from the large longwave CRF at the surface derived previously. Since negative values are 
possible but unlikely, one may conclude that the unreasonable high surface CRF is a result of the 
modeling and computation errors but not physical processes. Generally, the influence of clouds 
in the stratus regime causes the net longwave flux to be reduced by 50 - 60 W throughout 
the year. 

4. Discussion 

Monthly mean longwave radiation fluxes at the surface for four months have been 
determined from currently available satellite data. Because of the diurnal variation of surface 
temperature over the continents, surface net longwave fluxes, which involve the surface 
temperatures, are presented only over the oceanic areas. The method discussed here avoids the 
use of an independent estimate of the frequency of occurrence of clouds or even cloud top 
heights in determining the surface longwave fluxes. Current methods of modeling the longwave 
radiation processes in the atmosphere require certain assumptions regarding the presence of 
clouds and their horizontal and vertical extent. Because of the complicated nature of cloud 
radiation and uncertainties in the observadon of clouds, large errors are inevitable in this 
procedure. The procedure used in this study provides an alternative means of obtaining 
longwave radiative fluxes at the surface without the knowledge of cloud distribution. All 
meteorological parameters required in this method can be obtained from currendy available 
satellite data sets. The ISCCP data and U.S. Standard Atmosphere (COESA, 1976) provide the 
profiles ot temperature, water vapor mixing ratio and ozone as well as sea surface temperatures. 
The distributions ot the longwave CRF at the top ot the atmosphere are retrieved from ERBE 
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data which is currently being completed for several years of measurements. 

The weakest link in this procedure is the use of simulated CRF ratios that relate the CRF 
at the top and surface. Errors in the cloud generation scheme of the model used will atfect the 
ratio and hence, the final product. It is also not feasible to verify these ratios observationally on 
a global scale, although it would be useful to verify the model result in some specitic regions 
where simultaneous observations at the top of the atmosphere and the surface are available over 
an extended period of time. 

The final product can be compared with field data for a few regions and periods during 
which extended time observations were made from ships. Note that the fields generated are area 
averaged monthly means and comparisons with point measurements over short- periods of time 
do not provide any corroboration. Observations reported by Reed and Halpem (1975) off the 
Oregon coast covered an eight week period in July - August 1973 and included measurements 
taken from two sites, one 13 km from the shore and the other 120 km away. Lind and Katsaros 
(1987) have reported measurements taken off the California coast during the first two weeks ot 
November 1984. Radiation budget data from GATE (Cox and Griffith, 1979) was an amalgam 
of direct radiation measurements and modeled fields based on soundings. The authors point out 
the futility of computing area and time mean quantities directly from measurements even for a 
dedicated field campaign over a (relatively) small portion of the oceanic area of the globe. 

The daily mean net longwave radiative flux at the surface reported by Reed and Halpem 
(1975) varied from 71 W m'^ for days when the daily mean cloud cover was 10% - 20% to 
between 1 1 - 15 W m*^ for 70% - 100% cloud cover. An unweighted mean of the ten daily mean 
values taken at two stations over the period July 5 - August 26 is 33 W m The corresponding 
estimate for July 1985 from Fig. 8 is around 40 W m'^. This suggests that the procedure used 
here is able to provide a good measure of the mean cloud fraction which is the strongest 
determinant of the net longwave radiation tor this region. This bears out our thesis that a reliable 
measure of the effect of clouds is a necessary condition for obtaining global fields of the net 
longwave radiation at the surface. 

Lind and Katsaros (1987) have reported ship based observations of the upward and 
downward longwave radiation taken off the California coast from October 30 - November 14 
1984. The daily mean net longwave radiation from R/P FLIP ranged from 11 W m'2 to 69 W 
m*^. Although there is no reported cloud cover, inspection of the daily mean insolation shows 
that the extremes correspond to days of complete cloud cover and essentially clear skies, 
respectively. The mean for the 15 day period at the single station was 45 W m'^. The mean for 
October 1985 and January 1986 from Fig. 8 is 75 W m"^- This compares favorably with the 
satellite and model based results of Wu and Cheng (1989) for January 1979. The discrepancy 
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with the point measurement is significant but it is difficult to draw any conclusions based on the 
information available. From Fig. 3, the clear sky downward longwave flux for this region is 
about 315 W m‘^ in October 1985 and 275 W m‘^ in January 1986. The data from Lind and 
Katsaros (1987) indicate that the minimum daily mean value over the period was 340 W m'“, 
presumably on the clearest day. Inspection of the temporal data shows that the absolute 
minimum value reached was 300 W m'- on the night of November 7-8. This indicates that the 
clear sky element of the procedure is acceptable. The discrepancy is in the cloud estimate and 
there is no reasonable means of comparing the point measurement with an area average (for a 
different year). This highlights the difficulties inherent in validating global fields for a quantity 
that is measured only occasionally. Wu and Cheng (1989) have published maps of the 
downward and net upward longwave fluxes for January and July 1979 obtained by their 
physically based satellite retrieval. Their maps for July 1979 are reproduced in Fig 9. 
Comparison with the July 1985 panels in Fig. 4 and Fig. 8 show similarities and differences but 
the main features are common and can be considered to be the July climatology of the longwave 
fluxes at the surface. Wu and Chang (1991) have compared retrieved moisture, temperature and 
cloud cover from five sources and conclude that differences in the downward flux at the surface 
of 15 - 30 W m'2 should be expected when computations are made based on data from these 
different sources. 

Some further generalizations can be made about the uncertainty in the global fields. 
Fung et al. (1984) state that a 1 g kg'^ uncertainty in the water vapor mixing ratio at low levels 
and a 100 mb error in the cloud base each translate to a 10 W m'^ uncertainty in the net 
longwave flux over oceans. The presence of low level haze which is only now being mapped 
from space based measurements (Rao et al., 1989) and is not considered in the GCM is 
equivalent to underestimating low cloud cover, hence surface forcing. In certain regions this 
error could range from 5 - 10 W m'^. Moreover, if one makes the gross assumption that low 
clouds cover the oceans 50% of the time everywhere, it may be shown that the net longwave at 
the surface will range from 40 - 60 W m‘^ for any reasonable temperature and mixing ratio 
profile. The departure from this range of values shown in Fig. 8 is the true information 
contained in the maps and reflects the satellite based cloud information that has been used. 

5. Conclusions 

An attempt has been made to produce monthly mean global fields of the net longwave 
radiation flux at the surface over the oceans without resoiting to direct measurements. The key 
ingredients in the technique are the satellite derived temperature and moisture profiles (which are 
available operationally), the top of the atmosphere cloud radiative forcing from ERBE (which is 
an experimental data product but could be available in the future) and the cloud distributions 
from a general circulation model. It should be stres.sed that the actual cloud cover generated by 



13 


the model (which is subject to a great deal of uncertainty) is not used directly, hut only 
information on cloud type is used through a ratio of the cloud radiative forcing at the top and 
surface. This parameter has the advantage that it involves simulated radiation tluxes and not the 
cloud fraction. The latter quantity may vary considerably from model to model and also be quite 
different from satellite derived estimates. However, the radiation parameterization is usually 
such that the fluxes at the top of the atmosphere simulated by these models compare quite 
favorably with observations. Examples of this apparent contradiction are in Harshvardhan et al. 
(1989), and Kiehl and Ramanathan (1990). 

Although it is preferable to map global fields of the longwave radiation at the surface 
using space and ground based measurements alone, it is evident that the introduction of large 
scale numerical models into this effort is unavoidable. The surface is inaccessible to space based 
instruments at these wavelengths under cloudy conditions. Moreover, since the fluxes are 
extremely sensitive to water vapor mixing ratios near the surface, even clear sky estimates are 
subject to a great deal of uncertainty. Any global surface measurement program over the oceans 
is impractical. Hybrid techniques such as the one reported here that use several sources of data, 
both real and simulated, are the only options. 
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Figure Captions 

Flow diagram of the procedure to obtain maps of the monthly mean net upward 
longwave radiation flux at the surface using infonnation provided by a general 
circulation model (GCM), data from the Earth Radiation Budget Experiment 
(ERBE) and the International Satellite Cloud Climatology Project (ISCCP). 

Monthly mean longwave cloud radiative forcing at the surface for April, July, 
October 1985 and January 1986 obtained from ERBE top-of-ihe-atmosphere 
cloud forcing and GCM simulations of cloudiness. 

Monthly mean clear sky downward longwave radiation flux at the surface for 
April, July, October 1985 and January 1986 obtained from TOYS profiles on the 
ISCCP C 1 tapes and a broad band radiation code. 

Monthly mean atmospheric downward longwave radiation flux at the surface tor 
April, July, October 1985 and January 1986 obtained from the clear sky values 
shown in Figure 3 and the cloud forcing shown in Figure 2. 

Zonal mean sea surface temperatures (SST) for January 1986 from three different 
parameters on the ISCCP C 1 tapes and the blended SST from the Climate 
Analysis Center (CAC). 

Difference between the monthly mean surface emission in W m'- for January 
1986 computed using the CAC SST and the VIS/IR clear TS on the ISCCP Cl 
tapes. A positive difference indicates that the surface emission implied by the 
CAC SST is higher. 

Monthly mean clear sky net upward longwave radiation flux at the ocean surface 
for April, July, October 1985 and January 1986. 

Monthly mean atmospheric net upward longwave radiation flux at the ocean 
surface for April, July, October 1985 and January 1986. 

Monthly mean atmospheric downward longwave radiation flux at the surface (a) 
and the net upward longwave radiation flux at the surface (b) for July 1979 from 
Wu and Cheng (1989). 
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ABSTRACT 

Computational results have been obtained for the spherical albedo, global transmission, and global absorption 
of plane-parallel layers composed of cloud droplets. These computations, obtained using the doubling method 
for the entire range of single scattering albedos (0 ^ <i>o < 1 ) and for optical depths between 0.1 and 100, are 
compared with corresponding results obtained using selected multiple scattering approximations. Both the 
relative and absolute accuracies of asymptotic theory for thick layers, three diffuse two-stream approximations, 
and two integrated two-stream approximations arc presented as a function of optical thickness and single scanering 
albedo for a scattering phase function representative of cloud droplets at visible wavelengths. The spherical 
albedo and global absorption computed using asymptotic theory are found to be accurate to better than 5% for 
all values of the single scanering albedo, provided the optical thickness exceeds about 2. The diffuse two-stream 
approximations have relative accuracies that are much worse than 5% for the spherical albedo over most of the 
parameter space, yet arc accurate to within 5% in the global absorption when the absorption is significant. The 
integrated delta-Eddington scheme appears to be the most suitable model over the entire range of variables, 
generally producing relative errors of less than 5% in both the spherical albedo and global absorption. 


1. Introduction 

The role of clouds in determining the earth’s radia- 
tion budget has led to increased interest in the param- 
eterization of the radiative properties of cloud layers 
in numerical atmospheric models. Recent work has 
been concerned with relating cloud microphysics to 
optical properties (Slingo 1989) that can then be used 
in radiative transfer schemes within models. Most 
models now use some form of approximation to com- 
pute cloud radiative properties, such as the plane albedo 
from a given set of optical properties (optical thickness, 
sin^e scattering albedo, etc.). Whereas in the past these 
optical properties were generally fixed, there is now 
increasing use of interaaive schemes in which cloud 
optical properties are generated internally by the model 
(Charlock and Ramanathan 1985; Harshvardhan et 
al. 1989). 

As cloud fields evolve during a model integration, 
the optical properties of the generated clouds and 
models of gaseous absorption are used in a radiative- 
transfer scheme to provide the shortwave and longwave 
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radiative-energy field through the atmosphere. These 
computations need to be carried out at each model 
grid point at least every time the model cloud fields 
are updated. In models that resolve the diurnal cycle, 
this could be every three hours of simulated time, or 
even hourly. The computational burden is such that 
rapid, yet accurate, techniques are essential. In the 
shortwave, a common procedure is the computation 
of cloud-layer properties by a two-stream method and 
the adding of radiative fluxes through the atmosphere 
in an energy-conserving scheme (Lacis and Hansen 
1974; Coakley et al. 1983; Charlock and Ramanathan 
1985; Geleyn and Hollingsworth 1979; Harshvardhan 
et.al. 1987), although the two-stream equations can 
also be solved directly for multiple layers using matrix 
solvers (Wiscombe 1977; Toon et al. 1989), The flux 
adding method is essentially a severely truncated form 
of the adding-doubling method (Hansen and Travis 
1974), using upward and downward fluxes instead of 
intensities. 

In order to compute radiative fluxes through several 
atmospheric layers by the flux adding method, the ra- 
diative properties of cloud layers for two different 
sources are required (Harshvardhan et al. 1987; Kiehl 
et al. 1 987 ) . When collimated solar radiation is incident 
on an isolated cloud layer at some zenith angle with 
respect to the vertical direction, the fluxes emergent 
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from the layer in the upward and downward directions 
are determined by the plane albedo and total trans- 
mission of the layer, respectively. If the incident source 
is diffuse, the emergent flux may be obtained by an 
angular integration over the incident intensity field. In 
two-stream methods, the angular distribution of the 
incident intensity field is not resolved, and a common 
practice is to assume an isotropic diffuse source. For 
example, in a multilayer cloud system, the diffuse solar 
flux transmitted through the upper layer is the incident 
source for the lower layer. Also, in the case of a cloud 
layer overlying a reflecting ground surface, multiple 
reflections between the cloud and ground are consid- 
ered by assuming an isotropic diffuse source at the bot- 
tom boundary of the cloud layer. These diffuse radiative 
properties have also been used in the past to provide 
estimates of global effects of aerosol layers (Chylek and 
CoakJey 1974). A comprehensive study of the accuracy 
of various multiple scattering approximations for the 
plane albedo, total transmission, and fractional ab- 
sorption of isolated cloud layers corresponding to in- 
cident collimated radiation was presented by King and 
Harsh vardhan ( I986a,b). The present study comple- 
ments the earlier one in assessing the accuracy of var- 
ious approximations for calculating the radiative prop- 
erties of cloud and aerosol layers for an incident iso- 
tropic diffuse source. 

The presentation follows the organization of King 
and Harsh vardhan ( 1 986a, hereafter referred to as 
KJi). Section 2 discusses multiple scattering calcula- 
tions used to obtain the diffuse radiative properties of 
cloud layers of varying optical thicknesses and single 
scattering albedos. These computational results, ob- 
tained with the doubling method, will be considered 
the benchmark solutions with which various multiple 
scattering approximations will be compared. Section 
3 introduces the asymptotic theory approximation and 
the general class of two-stream approximations that we 
will consider. Section 4 presents the results of the com- 
parison between the approximate and exact results in 
terms of absolute and relative differences. A discussion 
of the results follows in section 5. Section 6 is a sum- 
mary including recommendations for using these ap- 
proximations. 

2. Multiple scattering computations 

To assess the accuracy of various multiple scattering 
approximations, radiative transfer computations were 
performed using the doubling method described by 
Hansen and Travis (1974), together with the invariant 
embedding initialization described by King (1983). 
These computations were performed for a cloud drop 
size distribution typical of fair weather cumulus ( FWC) 
clouds ( Hansen 1971), and were performed at a wave- 
length X = 0.754 /xm assuming a refractive index of 
liquid water m = 1,332. A detailed description of the 
cloud model, together with an illustration of the single 
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scattering phase function, can be found in KH. The 
azimuth-independent terms of the reflection and 
transmission functions were used to obtain the plane 
albedo r(r^, plq) and total transmission r(r^, as a 
function of the total optical thickness of the layer, 
and fjLo, the cosine of the solar zenith angle. In terms 
of these functions the spherical albedo, global trans- 
mission, and global absorption of the laver are given 
by 


r{T() - 2 1 r(r„ 

(1) 

Jo 

7(rr) = 2 I r(r,, 

Jo 

(2) 

d(T,)= 1 - F(t,) - 7(r,);- 

(3) 


In order to cover a wide range of applications, these 
computations were performed for values of the single 
scattering albedo ranging from pure absorption (a;o 
= 0) to conservative scattering (cjq = I ). The single 
scattering phase function was left unchanged such that 
all computations apply to a phase function having an 
asymmetry factor g = 0.843. 

Figure 1 illustrates numerical computations of the 
spherical albedo [r(rj, global transmission [7 (t;)], and 
global absorption [d(rr)] as a function of wq and r^. 
The doubling computations used to generate these re- 
sults were obtained at 12 optical depths 0,0625, 0. 125, 

. ♦ . , 128 interleaved with another set of 1 1 optical 
depths 0.0884, 0. 1 768, . . . , 90.5 1 . Each set of doubling 
computations was itself made at each of 31 values of 
the single scattering albedo. The single scattering albedo 
scale is linear in the similarity parameter s, defined by 



This makes it possible to expand the scale in the vicinity 
of conservative scattering ( o>o = I ) and still to span the 
full range 0 < wq ^ 1. The angular computations, in- 
cluding the integration in ( I ) and ( 2 ) , were performed 
at 80 Gaussian quadrature points. As in KH, the com- 
puted results were first interpolated to generate a 300 
X 300 matrix prior to plotting. The interpolated arrays 
represent the exact results to which the radiative trans- 
fer approximations are compared in section 4. 

It is perhaps pertinent to point out certain features 
of the radiative properties illustrated in Fig. 1 . For con- 
servative or very weakly absorbing layers, the spherical 
albedo increases rapidly with increasing optical thick- 
ness for small values of r, and then much more slowly 
as T( becomes large. This is the well-known nonlinear 
behavior that leads to problems in estimating area-av- 
eraged albedos for a nonhomogeneous cloud layer 
( Harsh vardhan and Randall 1985). For moderate to 
strong absorption, the saturation of both the spherical 
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a) 



c) 



albedo and global abso^tion at optical thicknesses of 
about 10 or even less is the most striking feature of 
Fig. 1. In the near-infrared, this implies that cloud ab- 
sorption is primarily a function of the single scattering 
albedo and not the optical thickness once the cloud 
layer is several hundred meters thick (Twomey 1976). 
The imponance of determining the spectral depen- 
dence of cvo for cloud layers and the development of 
accurate parameterizations for inclusion in radiative 
transfer models follows from this observation ( King et 
ai. 1990; Fouquan et ai. 1991). 



Optical Thickness 


Fig. 1. Doubling computations of the (a) spherical albedo 
Wo), (b) global transmission, 7(t,, wo), and (c) global absorption 
a(r,, t**o) as a funaion of optical thickness and singlc-scancring 
albedo for the FWC phase funaion. The single scattering albedo 
scale is linear in the similarity parameter, defined by Eq. (4). 


3. Radiative- transfer approximations 
Three classes of approximations will be considered 
here for comparison with the multiple scattering results 
presented above. In all cases, analytic or easily integra- 
ble functions relate the radiative properties to the 
optical properties. The three approximations we will 
consider are asymptotic theory for thick layers, diffuse 
two-stream approximations, and integrated two-stream 
approximations. Although there are several variations 
of two-stream approximations, only a few common 
and representative models will be considered. 
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a. Asymptotic theory 

Asymptotic theory is a rigorous solution to the 
equation of transfer in optically thick layers, and as 
such, makes no assumption about the angular distri- 
bution of scattered radiation within the medium. 
Expressions for the plane albedo and total transmission 
of an optically thick layer under collimated illumina- 
tion conditions can be found in KH and will not be 
repeated here. From these expressions, it can be shown 
that the asymptotic theory approximations for the 
spherical albedo [r(r^)], global transmission 
and global absorption (a(r^)] are given by 


set of differential equations for the upward and down- 
ward diffuse fluxes F-{t) (Coakiey and Chylek 1975; 
Meador and Weaver 1980) 

dF-{r) 

— — — = y\F {r)-yiF (r), (10) 

dF^ir) 

—^ = y,F-{T)-y,F^{r), (II) 

where F~{r) represents the upward flux and F'^(t) the 
downward flux at optical depth r. The equations can 
easily be solved subject to the boundary conditions 


r(rt) = 


1 — ’ 


(5) 


F^{0) = Fo, (12) 

F“(r,) = 0, .. (13) 


^ mn^e ' 

1 _ llg-2kr, ’ 

(6) 

a(T,) = 1 - r(r,) - 7 (t,), 

(7) 


for nonconservative scattering (ojq < I). In these 
expressions, is the spherical alb^o of a semi-infinite 
atmosphere and m, /t, /, and k are constants (coeffi- 
cients) that depend primarily on the similarity param- 
eter given by (4). All of the functions and constants 
that appear in these expressions can be computed by 
equating asymptotic formulas and doubling results at 
three values of the optical thickness for which asymp- 
totic theory is valid (viz., r, = 8, 16, and 32), as first 
pointed out by van de Hulst ( 1 968 ). Similarity relations 
for calculating (denoted ^ ^ by van de Hulst 1968 ), 
rt, /, and kzs 2 l function of s for the full range 0 
< j ^ 1 can be found in Table I of King et al. ( 1990). 
Once these coefficients have been computed, expres- 
sions for all of the radiative properties are analytic 
functions that can be computed rapidly within a ra- 
diative transfer code. 

For the special case of conservative scattering (cjq 
= 1 ), Eqs. (5) and (6) reduce to 


Hrt) = I 


4 

3(1 -^)(r, + 2^o) ’ 


( 8 ) 




4 

3(1 -^)(r, + 2^o) ’ 


(9) 


where Qo is the extrapolation length. The reduced ex- 
trapolation length ^' = ( I — ^)^o is known to range 
between 0.709 and 0.715 for all possible phase func- 
tions (van de Hulst 1980), and has the value q' = 0.715 
for the phase function used here. Again, one is left with 
simple ai\alytic functions describing the variation of 
r(r,) and /(r^) as a function of for a given asymmetry 
factor g. The set of equations (5)-(9) forms the ap- 
proximations for the diffuse radiative properties of a 
medium based on asymptotic theory. 


b. Diffuse two-stream approximations 

In the absence of any direct collimated beam, the 
two-stream equations of radiative transfer result in a 


for a diffuse isotropic source incident at the top bound- 
ary of the layer (or cloud) and for which no illumi- 
nation is incident from below. The spherical albedo is 
thus obtained from the expression 

Krt) = F’(0)/Fo, (14) 

and the global transmission from 

rir,) = F"(r,)/Fo. (15) 

For nonconservative scattering (ojo < 1 ), the solution 
may be obtained in the form (Coakiey and Chylek 
1975; Meador and Weaver 1980) 


= 72( 1 - g ^") 

^ + 7i + (* - 


(16) 


and for conservative scattering (cjq = 1 ) 

= (18) 
1 + y\rt 

/(r,)= 1 -r(r,). (19) 

In (16)-(19), the coefficients 7 1 and 72 depend on 
the particular two-stream approximation, with the dif- 
fusion exponent k defined as 


^ = (7i - 72 )'^^ 


( 20 ) 


Table 1 lists three diffuse two-stream models used 
for this study and the corresponding values of 71 and 
72* The discrete ordinates model is identified as the 
quadrature scheme by Meador and Weaver ( 1980) and 
Toon et al. ( 1989). The hemispheric-mean model de- 
fined by Toon et al. ( 1989) is similar to the Coakley- 
Chylek model II referred to by KH and first introduced 
by Chylek and Coakiey ( 1974). The two-stream mode! 
used by Sagan and Pollack ( 1967) has coefficients sim- 
ilar to those of both of the aforementioned models. 
Instead of the asymmetry parameter g, some two- 
stream models use the average backscatter fraction 
which is defined in KH and readily computed from 
the backscatter fraction /5(/i^), introduced by Coakiey 
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and Chylek ( 1975) and Zdunkowski et aJ. ( 1980) to 
compute the radiative properties of layers for colli- 
mated incident sources. The Eddington model has, of 
course, been used widely (Shettle and Weinman 1970). 
The set of equations ( l6)-(20) is used to compute the 
diffuse radiative properties for the two-stream approx- 
imations. It should be noted that these expressions have 
fairly simple analytic forms that favor rapid compu- 
tation. 

c. Integrated two-stream approximations 

Extensive discussion of two-stream approximations 
for a collimated source can be found in KH as well as 
in earlier work, in particular the comprehensive treat- 
ment by Meador and Weaver ( 1980). Expressions for 
the approximate plane albedo [r(r,, /xo)], total trans- 
mission [?(r^, Mo)], and fractional absorption [a(r,, 
Mo)] are the set of equations (21 )-(29) in KH. These 
expressions include the transformations that are re- 
quired in the case of delta scaling (Joseph et al. 1976). 
To obtain comparable expressions for the diffuse ra- 
diative properties, r(r,, mo) and r(r^, mo) must be in- 
tegrated according to Eqs. ( 1 ) and ( 2 ) . These expres- 
sions, however, are quite complicated, and thus inte- 
gration in a closed form is not generally practical. 

An analytic expression for the spherical albedo in 
the Eddington and delta-Eddington approximations 
has been obtained by Wiscombe and Warren ( 1980) 
and involves exponential integrals that are not con- 
ducive to rapid computation within a model. For this 
study, r{rt, mo) and r(r^, mo) obtained by the delta- 
Eddington approximation were numerically integrated 
to provide F(t/) and 7(r,). King and Harsh vardhan 
found that the delta-Eddington approximation for col- 
limated illumination conditions is quite accurate over 
a wide range of r, and mo, especially when a;o is near 
unity. A model that performs well for optically thin 
layers over the limited range of o;o studied by KH is 
the plane albedo scheme of Coakley and Chylek 
( 1975), designated Coakley-Chylek model I by KH. 
Two-stream methods for collimated sources require a 
third coefficient, 73 , which appears with the source term 
and is thus not included in Eqs. (10) and (II). The 


expressions for 73 used by the two integrated models 
presented here are given in Table 1. 

The integrations in Eqs. ( I ) and (2) required to ob- 
tain the diffuse properties are performed using 80-point 
Gaussian quadrature, and the results should be con- 
sidered identical to an analytic solution for all practical 
purposes. The general form of the quadrature sum- 
mation is 

N 

r{T,) = 2 Y. ( 21 ) 

(~\ 

where m/ are the Gaussian quadrature points on the 
half space and w, are the corresponding Gaussian 
weights. This detailed integration, however, is of no 
practical value because the computational burden is 
onerous when applied to a global climate model. We 
have, therefore, also included results for the delta-Ed- 
dington and Coakley-Chylek ( I) models integrated us- 
ing two-point and four-point quadrature, respectively. 
The diffuse radiative properties can then be obtained 
with a computational effort comparable to that re- 
quired to compute properties for collimated radiation. 

4. Results 

We have examined both the absolute and relative 
accuracies of the spherical albedo, global transmission, 
and global absorption as a function of r, and ojo for the 
asymptotic approximation, as well as the Eddington, 
discrete ordinates, and hemispheric-mean diffuse two- 
stream approximations. Other diffuse two-stream ap- 
proximations that we have examined generally yield 
somewhat poorer results when compared to our dou- 
bling benchmark calculations. In addition, we have 
considered the integrated delta-Eddington and Coak- 
ley-Chylek ( I ) approximations computed using both 
80 points and a limited number of Gaussian quadrature 
points for integration over the solar zenith angle. 

Figure 2 illustrates a 4 X 3 plot composite of results 
for the absolute difference in the spherical albedo, 
global transmission, and global absorption for four of 
these models, where the first row applies to asymptotic 
theory and succeeding rows to the Eddington, discrete 


Table 1. Summary of 7 , coefficients in selected two-stream approximations. 


Method 

7i 

72 

73 

Diffuse 

Eddington 

1/4(7 - u;o(4 + 3^)] 

-1/4(1 - «o(4 - 3j)l 

— 

Discrete ordinates 

V3/2[2 1 + ^)| 

1 - g)\ 



Hemispheric mean 

2 - 1 +1) 


— 

Integraied 

Delta-Eddington 

1/4(7 - wy4 + 3j')l 

-1/4(1 -wi(4-3g')J 

1/4/2 - 3g'^) 

CoakJey-Chyiek (I) 

(1 ~ (iJo[ I ■“ i3(Mo)l|/Mo 




* { I - 5”' W( i 
g‘ =<?/(! ^ g) 




Single Scattering Albedo Single Scattering Albedo Single Scattering Albedo Single Scattering Albedo 
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ordinates, and hemispheric-mean approximations. In- 
dividual plots in the first column of Fig. I represent 
absolute errors in the spherical albedo, defined as 

Af(T|, ciJo) = r{Tt, Wo) ““ r{Ti, wq), (22) 

with succeeding columns representing corresponding 
errors in global transmission [a7(t;, wq)] and global 
absorption [AJ(r,, wq)]. With these definitions, posi- 
tive (negative) errors indicate that the radiative transfer 
approximation overestimates ( underestimates) the ex- 
act solution, taken as the computational results pre- 
sented in Fig. 1. The relative errors in the spherical 
albedo, global transmission, and global absorption are 
presented in Fig. 3, and are given in percent. It is nec- 
essary to consider the performance of a particular 
model in both a relative and an absolute sense in order 
to delineate a range of acceptability. 

Individual contour plots in Figs. 2 and 3 have been 
shaded to draw attention to those regions of greatest 
accuracy. For example, asymptotic theory is seen to be 
accurate to within 5% in reflection and absorption for 
T, ^ 2 and for all values of wq. In transmission, relative 
errors exceed 5 % for wq < 0.90 and 2 < r, < 8, but the 
absolute errors are so small (<0.03) that the approxi- 
mation could probably still be used without serious 
adverse results. It is evident from Figs. 2 and 3 that the 
asymptotic approximation provides accurate results for 
all three diffuse radiative properties over the entire 
range of wo as long as r, ^ 2. 

The three diffuse two-stream models considered here 
are seen to yield unacceptable errors in one or more 
of the radiative properties over regions that would nor- 
mally be encountered in modeling applications. Al- 
though the range of acceptability will depend on the 
particular application, one can consider a 5 % error in 
the spherical albedo as a standard for comparison. The 
spherical albedo is usually the parameter of choice in 
estimating the sensitivity of any radiative perturbation. 
When the value itself is small, however, an absolute 
error criterion is more useful. For optically thin layers, 
the absolute errors in spherical albedo are generally 
less than 0.0 1 for the discrete ordinates and hemi- 
spheric-mean approximations. Errors in global trans- 
mission are similar for all three models, while the Ed- 
dington and hemispheric-mean models are successful 
in estimating the global absorption of a layer when wq 
^ 0.99 and ^ 10 with errors of less than 1%. If the 
range of acceptability is relaxed to 5%, then the Ed- 
dington and hemispheric-mean models can be used 
for absorption when wq is as low as 0.95 except for 
optically thick layers. This covers the range of single 
scattering albedo encountered in water clouds 
throughout the visible and near-infrared spectrum 
(King et al. 1990). 


The two integrated two-stream methods studied in 
this investigation provide more accurate results for all 
three diffuse radiative properties as shown in Figs. 4 
and 5. The delta-Eddington model was shown by KH 
to be highly successful in estimating the plane albedo 
for conservative scattering. There was a marked deg- 
radation of performance when nonconservative cases 
were considered. The present study shows that this 
model, when integrated over an isotropic diffuse in- 
cident source, provides excellent results tor the spher- 
ical albedo and global transmission over most of the 
range of r, and wo- Errors in excess of 10% in global 
absorption are present for moderate optical depths (0.5 
< T/ ^ 5) when wq exceeds about 0.95. It may be seen 
from Fig. 4, however, that the absolute errors in global 
absorption are less than 0.02 throughout this region. 
In addition, the large relative error in global transmis- 
sion for optically thick absorbing layers is irrelevant 
since the global transmission is itself close to zero, as 
is the absolute error. The Coakley-Chylek ( I ) model 
provides results of comparable accuracy for optically 
thin layers. This is not surprising since KH showed 
that it was the most accurate of the two-stream models 
for this case. The delta-Eddington model, however, 
when integrated over all incident angles, is nearly as 
accurate as the Coakley—Chylek ( I) model for optically 
thin layers. Moreover, the accuracy of the integrated 
delta-Eddington model does not degrade as rapidly at 
higher optical depths. 

As mentioned previously, these two models would 
only be of academic interest if a rigorous numerical 
integration were required for every computation of the 
diffuse radiative properties. We have, therefore, also 
presented results obtained using a limited number of 
quadrature points in the integration over solar zenith 
angle [cf. Eq. (21 )] . As can be seen from the second 
panel of Figs. 4 and 5, a two-point integration of the 
delta-Eddington models yields accuracies that are 
comparable to the accuracy obtained using an 80-point 
integration. For the Coakley-Chylek (I) model, how- 
ever, it is necessary to use a four-point integration to 
obtain results that are of comparable accuracy. 


5. Discussion 

Although the results presented here are not exhaus- 
tive in the sense that all possible approximations have 
not been tested, we feel they are representative of what 
one might expect for any class of model. All the 
schemes are computationally efficient, and it is not 
necessary to perform a rigorous integration for the 
models based on incident collimated sources. The ap- 
proximations presented here can be incorporated into 


Fig. 2. Absolute accuracv ofasvmpiouc theory, Eddington, discrete ordinates and hemisphcnc-mcan approximations to the 
sphencal albedo, global transmission and global absorption as a function of optical thickness and single scaiienng albedo. The 
phase function is assumed throughout. 
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Fig. 5. As in Fig. 4 e.xcept for relative accuracies ( in percent). 
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a multilayer radiative transfer module that may be 
added to the radiation code in a numerical model. 

All computations presented here were obtained for 
a FWC drop size distribution having an asymmetry 
factor g = 0.843. Variations along the o>o axis can 
therefore be viewed as representing the etfect of altering 
the gaseous absorption in the layer at a particular 
wavelength, or to some extent, variations in the wave- 
length if ^ does not vary too greatly. This would cover 
the solar near-infrared spectrum over which (I ^ ojq) 
varies by several orders of magnitude, while ^generally 
lies between 0.80 and 0.90 (cf. King et al. 1990). 

As found by KH for collimated radiative properties, 
the asymptotic approximation yields consistently ex- 
cellent results for optically thick layers, regardless of 
single scattering albedo and solar zenith angle. Figures 
2 and 3 show that the same is true for the diffuse ra- 
diative properties as long as T; ^ 2. In a numerical 
model with internally generated cloud optical prop- 
erties, this requirement will not always be met. Errors 
become unacceptably large when < 1. For this rea- 
son, the asymptotic approximation should only be used 
when it is known a priori that r, ^ 2 at all times. This 
is the one serious shortcoming of an otherwise simple 
and accurate model. The method also requires a pre- 
computed table of coefficients m, n, k, /, and or 
analytic forms that compute these quantities within 
the program. Analytic expressions for these coefficients 
in terms of the similarity parameter can be found in 
King et al. ( 1990), which further discusses a remote 
sensing application of asymptotic theory. 

The three diffuse two-stream models presented here 
are the simplest to implement in a numerical atmo- 
spheric model and are the most computationally effi- 
cient, but their accuracy is limited to certain regions 
of the parameter space. They are also not uniformly 
accurate for all three radiative propenies. This is es- 
pecially true in the Eddington approximation, where 
the spherical albedo is frequently too inaccurate to be 
of any value in a numerical model. In addition, the 
Eddington model yields unphysical values of the 
spherical albedo and global absorption when absorption 
is very large (King and Harshvardhan 1986b). This 
situation arises occasionally in the water vapor bands 
of the near-infrared and frequently in the thermal in- 
frared. The problem can be rectified in a computer 
code with the addition of a check for unphysical values 
that could then be forced to the condition of zero re- 
flection. The discrete ordinates model does not suffer 
from this limitation and generally yields better results 
for the spherical albedo than does the Eddington ap- 
proximation, The somewhat poorer results for global 
absorption are not too important since the absolute 
errors are small in this case. The hemispheric-mean 
model yields results very similar to the discrete ordi- 
nates model, except for global absorption. The smaller 
relative errors for weak absorption are an especially 
attractive feature of the hemispheric-mean model. 


which otherwise suffers from the fact that it tends to 
overestimate the spherical albedo by more than 5% for 
the very important case of nearly conservative optically 
thick layers. 

The integrated delta-Eddington model yields excel- 
lent results for all three radiative properties over the 
entire range of optical properties that are encountered 
in the radiation code of a numerical atmospheric 
model. In fact, errors in the diffuse radiative properties 
are generally smaller than the errors found by KH for 
collimated radiative properties, with no unphysical re- 
sults anywhere in the parameter space. There has ob- 
viously been some cancellation of errors in the angular 
integration. As mentioned earlier, the one error-prone 
region is moderate optical thickness and weak absorp- 
tion. This was also true for the errors in fractional ab- 
sorption for a collimated source. Since the direct beam 
is usually handled by a delta-Eddington or similar ap- 
proximation, the coefficients and functions used for 
this model are usually already present in a numerical 
model. There is, however, an extra computational 
overhead in the angular integration, in that planar 
properties need to be computed at several angles and 
then numerically integrated. As seen in Figs. 4 and 5, 
however, these computations need be carried out at 
only two points to yield results comparable to a detailed 
numerical integration. 

The integrated Coakley-Chylek (I) model is of lim- 
ited value, except perhaps for optically thin, weakly 
absorbing layers. There is also an added computational 
burden since at least four angular computations are 
required for the phase function used here. For colli- 
mated radiative properties and for optically thin layers, 
KH found that this model was superior to the delta- 
Eddington model. For diffuse radiative properties, on 
the other hand, we find that there is little advantage in 
using the CoakJey-Chylek (I) model, even for optically 
thin layers. 

6. Summary and recommendations 

In the present study the spherical albedo, global 
transmission, and global absorption computed by var- 
ious radiative transfer approximations have been com- 
pared with doubling computations as a function of op- 
tical thickness and single scattering albedo. Since the 
entire range of ojo has been considered for optical depths 
from 0.1 to 100, the results presented here can be uti- 
lized to decide which approximate method is the most 
accurate for a particular application. The results pre- 
sented here should be considered in parallel with the 
findings of KH regarding the plane albedo, total trans- 
mission, and fractional absorption for a collimated in- 
cident source. 

In order to summarize the results of this study, it is 
useful to present composite figures extracted from the 
individual figures to highlight regions of highest ac- 
curacy. Following van de Hulst ( 1980), we show in 
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Fig. 6 the regions for which a particular model is ac- 
curate to within 1% and 5%. Only those models that 
are reasonably accurate in the particular radiative 
property have been included. These models include 
asymptotic theory, the two-point deita-Eddington 
method, and the four-point Coakiey-Chylek (I) 
method. Although the hemispheric-mean model yields 
acceptable results for the global absorption, it is not 
included here because results for the spherical albedo 
are generally poor. 

At the I % ( 5% ) level, asymptotic theory can be used 
for all (Vo as long as r, ^ 3.5 (2). For smaller optical 
depths, there is a choice that can be made between the 
deita-Eddington and Coakley-Chylek (I) models, but 
our recommendation is to use the deita-Eddington 
method. Many general circulation models are already 
using this method to compute collimated radiative 


properties, and the additional overhead incurred in the 
two-point integration should be minimal. If a scheme 
is needed to span the entire domain, the asymptotic 
method should not be used since its performance de- 
teriorates very rapidly for t, ^ 3. For this situation, 
typical of GCM applications, the integrated deita-Ed- 
dington scheme should yield acceptable results. 

The overall errors for a multilayer cloud system over 
a reflecting surface will depend on the optical thickness 
and single scattering albedo of the individual layers. 
At present, it is felt that errors in parameterizing the 
band-averaged single scattering albedo of cloud layers 
in the near-infrared will dominate errors in approxi- 
mating the radiative properties of individual layers 
( Fouquart et al. 1991). For example, the use of a single 
value of (Vo to represent the entire solar near infrared 
can result in errors in the layer absorption of several 


1% Error 




Optical Thickness Optical Thickness Optical Thickness 


Asymptotic Theory 

Two-Point integrated Deita-Eddington 

Four-Point Integrated Coakley-Chylek (I) 

Fig 6. T^c domain of validity of selected approximations to the diffuse radiative properties of cloud layers as a funaion of 
optical thicknca ( rj and single scattering albedo (wo). The upper panels correspond to a relative accuracy of 1% and the lower 
panels to 5,o. The single scattenng albedo scale is linear in the similarity parameter. The FWC phase function is assumed 
throughout. 
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hundred percent (Slingo 1989). The sensitivity of all 
radiative properties to ojq can be appreciated by in- 
spection of Fig. 1. Since any scheme has to limit the 
number of bands for computational efficiency, the se- 
lection of these bands and the average absorbing prop- 
erties used could determine the overall accuracy. For 
a given set of r, and a?o, however, the results presented 
in this study could act as a guide for choosing an ap- 
propriate model. Finally, it is pertinent to mention that 
these accuracies refer to an idealized plane parallel 
model. There is, of course, the additional problem of 
representing inhomogeneous cloud systems including 
geometric effects ( Harsh vardhan and Thomas 1984; 
Stephens 1988 ), a problem not considered in this study. 
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ABSTRACT 


Current atmospheric general circulation models parameterize cloud-radiation effects 
through relationships that are based on the assumpdon that a model grid is partly clear and the 
remaining area is filled with a homogeneous cloud layer. The optical properties of this homo- 
geneous layer are related to cloud microphysical quantities such as the liquid water path and 
effective drop radius which in turn are diagnosed from model variables. In the future it will be 
necessary to consider distributions of liquid water paths within a grid and relate the effective 
radius to the type of cloud formed since the grid mean radiative properties are very sensitive to 
these parameters. There is also the need for a computationally efficient technique to model 
overlapping liquid and vapor absorption in the solar near infrared. 


INTRODUCTION 

The role of clouds has been identified as the area of research of the highest priority in 
many global change programs (e.g. Committee on Earth Sciences. 1989). One reason is the 
dramatic radiative effect of clouds in both the solar and thermal regions of the electromagnetic 
spectrum. In general, clouds have opposite effects in these two spectral regions and the net 
effect for the earth-atmosphere column is the difference between these two large forcings. The 
global distribution of the individual and combined forcing at the top of the atmosphere has been 
studied extensively during the Earth Radiation Budget Experiment (ERBE). The ERBE results 
(Harrison et al., 1990) indicate that in the global annual mean, clouds have a cooling effect, i.e. 
the reflection of solar radiation dominates the thermal trapping effect. However, in areas such as 
the tropics, the two effects nearly cancel, but of course, the vertical profile of the forcing is quite 
different in the solar and in the thermal. 

In order to include the radiative effect of clouds in a numerical climate model it is es- 
sential to consider all the important processes that could change cloud cover and cloud radiative 
properties during a climate simulation. These are examples of interactive feedback processes 
that when included in climate models can alter their response to a given forcing such as doubling 
atmospheric carbon dioxide. Figure 1 shows schematically some possible feedback paths 
involving cloud microphysics. 

The most comprehensive tool for the simulation and analysis of climatic change is a gen- 
eral circulation model (GCM) of the atmosphere. Often, GCMs are coupled to simple ocean and 
sea ice models. The manner in which cloud-radiation effects are parameterized differs from 
model to model. This review will deal with the rationale behind these parametenzations and 

discuss possibilities for future improvements. 

RADIATIVE PROPERTIES OF HOMOGENEOUS CLOUDS 

All current GCMs have interactive cloud generation schemes and parametenzations for 
computing the radiative properties of cloudy layers. Earlier models fixed cloud radiative proper- 
ties (Wetherald and Manabe. 1980) which limited the ability of changes in cloud properties to 


3 


feed back to the climate system. An exhaustive list of cloud parameterizations is given in Cess 
et al. (1990) which reports on an intercomparison of 19 atmospheric GCMs that were subjected 
to a reverse climate experiment. The models were run twice in a perpetual July mode for 
prescribed sea surface temperatures (SST) that were uniformly 2 K above and 2 K below the 
climatological mean SSTs for July. All models showed decreased cloud cover for the warm 
simulation relative to the colder one but the climate sensitivity parameter (Cess et al., 1990), X,, 
whose magnitude depends on the various feedbacks in the models ranged over a factor of 3. 
This is because changes in the vertical distribution of clouds were different from model to model 
and the optical (hence, radiative) properties of clouds changed with the change in mean climate 
in a manner that was unique to each model. 

Cloud Microphysics 

The parameterization of cloud optical properties in terms of GCM resolved variables is a 
fairly recent trend and is based on first principles applied to homogeneous cloud layers. It may 
be shown that the extinction coefficient of a volume of cloud drops is related simply to various 
moments of cloud microphysical parameters in two limiting cases (Stephens, 1984). At solar 
wavelengths, the characteristic dimension of cloud drops is much larger (> 4 |im) than the wave- 
length and in this case the droplet extinction coefficient, a (km*^), is proportional to L/r^ where 
L (g m'^) is the liquid water content and r^ (jim) is the effective radius which is defined as the 
ratio of the third to the second moments of the drop size distribution (DSD). It is this inverse 
relationship to drop size that results in an increase in the extinction coefficient (and cloud 
albedo) when the DSD is shifted to smaller sizes (Twomey, 1977). 

The radiative properties of a homogeneous cloud layer are computed using the column 
integrated extinction coefficient which is a non-dimensional parameter called the optical depth. 
For example, the optical depth at visible wavelengths is 
top 

^vis ~ J*^vis ^ ^ ^ 

base 
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where z is geometric height. 

Since the column integrated liquid water path is 
top 

W(g m-2) = jX(g m-3)dz, (2) 

base 

then for visible wavelengths 

'Cvis W/rg^ (3) 

if r^ is uniform through the cloud. Furthermore, can be parameterized solely in terms of W 
if one assumes further that there is a systematic relationship between r^ and W (e.g. the model of 
Stephens, 1978). If can be related to W, cloud optical properties at other wavelengths 
follow immediately. For instance, the thermal flux emittance of a cloud layer is approximated 
by 

= 1 - exp (- (4) 

where cloud optical depth at thermal wavelengths is = 0.5 (Platt and Harshvardhan, 
1988). The diffusivity factor, p, which accounts for the angular distribution of thermal emission 
is roughly 1.5 so that the thermal flux emittance of a cloud layer may be written as 

e. = 1 - exp (- 0.75 x . ). (5) 

The above development can be extended to ice clouds with some modifications (Platt and 
Harshvardhan, 1988). 

Application to GCMs 

Equations (3) and (5) form the basis of most parameterizations of cloud optical properties 
in GCMs. The radiative properties of a cloud layer (reflectance, transmittance and absorptance) 
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are computed using a two-stream technique (King and Harshvardhan, 1986) which requires as 
input the droplet single scattering albedo, (O, and asymmetry parameter, g, of the scattering phase 
function in addition to X. Models generally prescribe these as wavelength dependent constants. 
The radiative properties of isolated cloud layers have to be integrated into the column model and 
this is usually accomplished through a flux adding (Harshvardhan et al„ 1987) or matrix 
inversion (Toon et al., 1989) technique. Application of the above procedure is possible only if a 
GCM provides information on the cloud water content and related microphysics. However, this 
aspect of modeling is in its infancy and many GCMs use the layer mean temperature as a 
surrogate for cloud microphysics (see Tables 4 and 5 of Cess et al., 1990). 

The validity of using temperature to parameterize cloud microphysics has not been 
thoroughly examined over the entire range of global conditions. However, there is empirical and 
theoretical evidence that optical properties can be related to temperature for cold clouds 
(Petukhov, et al., 1975; Feigelson, 1978; Heymsfield and Platt, 1984; Somerville and Remer, 
1984; Betts and Harshvardhan, 1987; Platt and Harshvardhan, 1988; Stephens et al., 1990). At 
the very least, the sharp change in effective radius accompanying phase change can be used to 
distinguish between water and ice clouds in a GCM (Charlock and Ramanathan, 1985). 
However, a direct relationship with temperature invariably breaks down if extended to 
convective anvils which are optically thick (Harshvardhan et al., 1989). The studies cited above 
also show that the optical properties of warm clouds cannot be parameterized in terms of 
temperature alone. In principle, it is possible to incorporate cloud formation and dissipation 
mechanisms into larger scale models (Heymsfield and Donner, 1990). 

MEAN PROPERTIES OF CLOUD FIELDS 

The horizontal grid increment of GCMs is of the order of 200 - 500 km within which one 
generally finds mesoscale structure as well as detailed structure at much smaller scales. There is 
need for subgrid scale parameterization of cloudiness in GCMs because of the stark contrast 
between clear and cloudy radiation fields as well as the non-linear behavior of radiative 
properties with respect to cloud optical properties. Most current GCM radiation routines 
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distinguish between clear and cloudy portions of the horizontal grid in each layer and they also 
include some scheme for overlapping fractional cloud cover in the vertical (Tian and Curry, 
1989). Tables 4 and 5 of Cess et al. (1990) list the cloud fraction generation schemes in several 
current GCMs. 

In principle, if there is some means of determining the cloud water (or ice) content in a 
grid volume, the application of the relationships given by eqs. (3) - (5) to a partially filled area 
can provide the overall grid box mean radiative properties of the layer. However, the 
distribution of the liquid water in the cloudy portion could affect the grid mean properties 
substantially. This was demonstrated by Harshvardhan and Randall (1985) for very simple, yet 
plausible, variations in the liquid water content of a spatial grid. Their results are reproduced in 
Figure 2 and show quite clearly that mean radiative properties of a cloudy layer can not be 
uniquely related to the mean liquid water column amount in the grid. 

There is also observational evidence of this non-unique dependence. Stephens and 
Greenwald (1991) have correlated the microwave derived liquid water path and albedo of clouds 
from an analysis of Nimbus 7 data. Figure 3 summarizes their study which shows that the 
albedo-liquid water relationships for homogeneous clouds derived from eq. (3) are applicable to 
clouds in midlatitudes but can not be applied in the tropics. In fact, the relationship in the 
tropics is much like Case C shown in Figure 2. There is obviously considerable structure in 
tropical clouds within the field of view analyzed in the study. In addition, there are cloud 
geometry effects which have not been considered in Figure 2. 

The current practice of separating a horizontal area into clear and cloudy fractions is a 
particular, though extreme, example of a subgrid scale liquid water distribution parameteriza- 
tion. Observational and modeling studies currently being conducted should lead to more general 
parameterizations. A possible source of information on subgrid scale variability is data from the 
International Satellite Cloud Climatology Project (ISCCP) which has cataloged cloud properties 
inferred from visible and infrared window radiances measured by instruments on geostationary 
and polar orbiting satellites (Rossow and Schiffer, 1991). The data coverage is nearly global 
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with a time resolution of three hours for the period July 1983 to the present. A parallel 
campaign of field observations under the FIRE (First ISCCP Regional Experiment) program 
(Cox et al., 1987; Albrecht et al., 1988) has also yielded a wealth of information on the 
horizontal variability of cloud microphysical properties in stratocumulus (Nakajima et al., 1991) 
and cirrus clouds (Spinhime and Hart, 1990). 

The three dimensional structure of cloud fields is currently being simulated by cloud 
scale models (Moeng, 1986; Tao et al., 1987; Xu and Krueger, 1991). Computational con- 
straints have prevented extended time integrations but in the near future, model realizations 
could be used to construct parameterizations for coarser resolution GCMs. The mean radiative 
properties of the cloudy portion of the field can then be computed from the distribution of 
model microphysical and macrophysical quantities and parameterized in terms of the mean and 
higher moments of, for example, liquid water and effective radius distributions. The model 
envisaged is shown schematically in Figure 4 for liquid water path. 

SPECTRAL RADIATIVE PROPERTIES 

The above discussion focused on computations of cloud radiative properties at visible 
wavelengths but of course is applicable for all radiative properties throughout the spectrum. At 
thermal wavelengths, the emittance of water clouds saturates at fairly modest thicknesses (see 
eq. 5) but a subgrid scale model as in Figure 4b would be useful for thin cirrus and even marine 
stratocumulus. One of the most pressing problems in cloud radiation modeling is the 
parameterization of cloud optical properties in the solar near infrared which contains almost fifty 
percent of the extra-terrestrial solar energy. At wavelengths greater than 0.7 |im, water in the 
vapor and condensed phases exhibits strong absorption features that overlap somewhat. Fur- 
thermore, cloud layer absoprtance depends critically on the single scattering albedo, co, which 
varies through the near infrared and, more importantly, depends on the DSD (King et al., 1990). 
The recent study of Intercomparison of Radiation Codes in Climate Models (ICRCCM) showed 
that low resolution models produced large errors because of the spectral averaging of co 
(Fouquart et al., 1991). 
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It is for the above reasons that solar radiation computations are made for several bands 
and the single-scattering properties are parameterized in terms of cloud microphysical quantities 
for each band separately. Slingo (1989) has proposed a four band model for which the optical 
properties, x, to, and g, are in the form 

Xi = (ai + bi/re) W (6) 

1 - cOj = Cj + djrg 


gi = ei + fjrg 

where aj - f. are band dependent coefficients that are obtained by fitting the parameterized 

results to more exact calculations and W and r^ are assumed to be given. Since models do not 

diagnose r^, the above relationships have been extended by Ritter and Geleyn (1992) to include a 

relationship of the form 

r = c + C-L 
e 1 2 

where c^ and c^ are coefficients. Ebert and Curry (1992) have obtained coefficients for eqs. (6) 
- (8) that are appropriate for ice clouds. Again, it is necessary to make some diagnosis of cloud 
microphysical properties in order to apply the relationships. 

The absorptive properties of condensed phase water in the near infrared cannot be con- 
sidered in isolation but must be integrated with water vapor absorption. This is usually accom- 
plished by modeling vapor absorption with the k-distribution technique (Lacis and Hansen, 
1974) which approximates the solution to the spectrally varying problem with an equivalent set 
of pseudo-monochromatic computations. Accuracy generally increases with the number of 
bands used and the number of k-values within bands. However, computation time rapidly 
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becomes prohibitive. Since liquid and vapor absorb at similar wavelengths, the spectrally 
integrated radiative properties of a cloud layer depend on the column amount of vapor above the 
cloud because the solar flux impinging on the layer does not contain energy that has already 
been depleted by the vapor above (Davies et al., 1984). Methods of incorporating this into a 
parameterization have been suggested by Fouquart and Bonnel (1980), Ramaswamy and 
Freidenreich (1992) and others, 

A suggestion by Briegleb (1992) that involves the identification of vapor absorption with 
different liquid absorption bands can reduce computation time significantly. Figure 5(a) shows 
the solar spectral column absorption by water vapor and a thick cloud of liquid droplets 
separately in the absence of the other constituent for a solar zenith angle of 30° and a standard 
midlatitude summer water vapor profile. The liquid DSD is from King et al. (1990) and has an 
effective radius of 8 (tm. The cloud is assumed to be thick enough such that radiative properties 
have attained their asymptotic values. It may be seen that there is a tendency for vapor and 
liquid absorption to increase in lockstep with increasing wavelength (decreasing wavenumber). 
This feature motivated Briegleb (1992) to assign the different k-values from Lacis and Hansen 
(1974) to distinct droplet absorption properties based on the Slingo (1989) parameterization. 

There is, however, a flaw in this argument that can be appreciated by inspection of 
Figure 5(b) which is a scatter plot of the two separate absorptions shown in Figure 5(a) sampled 
at 100 cm-1 intervals. Whereas it is true that high vapor absorption is associated with high liquid 
absorption, there is fairly strong liquid absorption even when vapor absorption is low. In fact, if 
there is a substantial column of vapor above the cloud, all the solar absorption in the cloud layer 
will occur in these window regions. The general idea of relating vapor and liquid absorption can 
still be applied by deriving a more quantitative relationship based on k-distributions at higher 
resolution. 
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An example of such a relationship made using the k-distribution given in Table 1 of 
Chou (1986) is shown in Figure 6. The k-distribution is solar flux- weighted and computed for 
p = 300 mb and T = 240 K. Cloud layer absorption is linearly proportional to the single 
scattering albedo, (D, in the thin limit and is primarily a function of the similarity parameter, s, 
where 

s = [(l -to)/(l -cog)]"'^, (10) 

for very thick clouds (King and Harshvardhan, 1986). By assigning the appropriate broad band 
liquid properties to each of the k-values in the particular band and computing the weighted sum, 
two separate relationships are obtained for the limiting cases of exceedingly thin and thick 
layers. The relationships shown in Figure 6 are unique to the cloud drop size distribution 
considered but a more flexible model could perhaps be a useful alternative to current techniques. 
If computation time is not a constraint, then a multi-band model is probably preferable. 
SUMMARY 

Since the early climate models of the 1960s there has been a gradual progression of 
increasingly interactive cloud radiation parameterizations. This has allowed the various cloud 
radiation feedbacks to operate in models. One result has been the wide range in climate sensi- 
tivity displayed by models having different parameterizations. The intercomparison study of 
Cess et al. (1990) has clearly demonstrated the divergence in the treatment of clouds in current 
models. 

Future thrusts will probably be in the areas of modeling and parameterizing cloud 
microstructure and also cloud field macrostructure. Climate models will continue to get more 
interactive and the role of clouds in model simulations will become increasingly important when 
oceanic coupling is included. Studies with cloud scale and mesoscale models are already 
influencing GCM parameterizations (Xu and Krueger, 1991). This trend will accelerate as the 
various groups cooperate in the study of the role of clouds in maintaining the current climate and 
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modifying the climatic response to anthropogenic forcings. 
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Figure Captions 


Fig. 1. Some plausible cloud feedback processes. 

Fig. 2. (a) Schematic illustration of three liquid water path (W, g m probability 

distributions considered as examples; (b) and (c) show the diffuse albedo and 
emittance, respectively, for the three cases illustrated in (a), after Harshvardhan and 
Randall (1985). 

Fig. 3. Ellipses containing 97% of the annual composite cloud albedo-liquid water path 
relationship data for the tropics and midlatitudes. Also shown are the results of 
theoretical calculations (solid curves) for plane-parallel clouds with various effective 
radii rg (lim), after Stephens and Greenwald (1991). 

Fig. 4. Schematic diagram of the liquid water path distribution in (a) current GCMs and (b) 
future models. 

Fig. 5. (a) Solar near infrared absorption by water vapor (thin line) and cloud droplets (thick 

line) separately in the absence of the other constituent for a midlatitude summer 
water vapor profile. The effective radius of the droplets is 8 ^irn and the cloud is 
assumed to be optically thick. Solar zenith angle is 30°. The outer envelope shown 
by the dashed line is the extraterrestrial solar flux. All fluxes have been averaged 
over 100 cm'^; (b) scatter plot of vapor and liquid absorption shown in (a). 

Fig. 6. Relationship between droplet single scattering albedo and vapor absorption for the 
limiting cases of thin and thick cloud layers. 
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